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DRSAR-SAM 


MEMORANDUM  FOR  RECORD 

SUBJECT:  User's  Guide  for  the  Computer  Program:  Exterior  Ballistics  of 
Boosted  Rockets  (EXBAL) 


1 . References : 

a.  Technical  Report,  DRSAR/SA/R-12,  April  76,  title:  Dynamics  of 
Liquid-Filled  Projectiles. 

b.  MFR,  DRSAR-SAM,  20  Apr  77,  subject:  An  Improved  Algorithm  for 
the  Glide  Mode  of  Copperhead  Used  in  a 3 DOF  Flight  Simulation. 

2.  A general  purpose,  point-mass  (3  DOF)  flight  simulation  program, 

EXBAL,  has  been  used  by  DRSAR-SA  and  others  for  many  years  for  a variety  of 
applications.-  A recent  special  application  is  noted  in  Ref  a.  Each  appli- 
cation has  special  requirements  which  often  entail  modifications  to  the 
basic  program.  Consequently,  the  program  documentation  must  be  updated 
frequently.  • 

3.  The  purpose  of  this  memorandum  is  to  update  EXBAL  for  the  program  user. 
Recent  changes  that  were  incorporated  to  treat  the  glide  mode  of  Copper- 
head (Ref  b)  are  discussed  here.  In  Annex  1 (Incl  1)  to  this  memo  are 
contained  the  following:  background,  program  structure,  theory  pertinent 
to  recent  changes,  input  and  output  data  definition  and  format,  flow  charts 
for  the  major  subprograms,  and  bibliographic  citations.  Annex  2 (Incl  2) 
contains  the  source  program  listing  in  Fortran  4 with  an  example. 


2 Incl 
as 


^ O--  <■  '•I 

GEORGE  SCHLENKER 
Operations  Research  Analyst 
Methodology  Division 
Systems  Analysis  Directorate 


Next  page  is  blank. 
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ANNEX  1 


DESCRIPTION  OF  THE  COMPUTER  PROGRAM: 
EXTERIOR  BALLISTICS  OF  BOOSTED  ROCKETS  (EXBAL) 


1.  General 

This  program  is  a multi-purpose,  point-mass  flight 
simulation.  It  was  conceived  as  a means  of  simulating 
trajectories  for  spin-  and  fin-stabilized  projectiles 
and  rockets  in  three-space.  The  principal  intended  ap- 
plication is  to  systems  which  can  be  assumed  aerodynami- 
cally  stable  so  that  trailing  or  following  behavior  of 
fin-stabilized  systems  is  exhibited  and  so  that  the  equi' 
librium  yaw  of  repose  for  spin-stabilized  systems  is 
quickly  achieved.  Important  considerations  in  develop- 
ing the  program  were: 

a.  operating  efficiency  or  speed  of  execution 

b.  ease  of  use  by  unsophisticated  users 

c.  minimal  input  data  requirements 

d.  ability  to  generate  multiple  trajectories  under 
program  control  for  parametric  analysis 

e.  ability  to  treat  a variety  of  system  types 
via  data  changes  and  option  switches 

f.  modularity  for  ease  of  program  modification. 

The  following  sections  will  treat  potential  program 
applications,  subprograms  employed,  execution  options, 
input  variables,  output  variables,  and  flow  charts. 
Documentation  for  most  of  the  theory  is  supplied  in  BRL 
Memorandum  Report  No.  I6l7,  September  1961j.  [2]  and  BRL 
Report  No.  131^»  March  1966  [li]. 

2.  Applications 

This  program  was  originally  developed  to  analyze 
conceptual  systems  such  as  fin-stabilized  gun-boosted 
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rockets  in  which  the  principal  interest  was  in  achievable 
range.  However,  subsequent  applications  involved  con- 
ventional, purely  ballistic  systems  in  which  significant 
variables  were  range,  maximum  ordinate,  deflection  due 
to  spin,  and  time  of  flight.  Other  applications  explored 
differential  effects  for  error  analysis  associated  with 
projectile  inertial  properties,  meteorological  condi- 
tions, and  launch  conditions.  Recently,  program  changes 
were  made  to  treat  projectiles  having  wings,  fins,  and 
controls  for  a glide  mode  of  flight. 

3.  Subprograms 

The  program  is  organized  modularly  into  subprograms 
each  of  which  performs  a separate,  discrete  fixnction. 

The  MAIN  program  performs  such  executive  functions 
as  input,  -control  of  plotting  and  printing  output,  change 
in  time  step  as  required,  stepping  thru  parameter  loops, 
and  termination  of  run. 

Subprogram  RUNGE  contains  a fourth-order  Runge-Kutta 
integration  algorithm.  The  order  of  the  equations  is 
specified  as  is  the  number  of  independent  variables  in 
an  entry  point  RUNGEl.  The  initial  values  of  higher 
derivatives  are  computed  on  the  occasion  of  the  first 
call.  Subsequent  subroutine  calls  are  made  to  entry 
point  RUNGE2  at  which  the  system  state  is  updated  from 
t to  t + dt. 

Subprogram  SOUND  contains  all  atmospheric  data. 

This  program  generates  the  local  speed  of  sound,  local 
gravitational  acceleration  (exclusive  of  centrifugal 
acceleration) , air  density  and  viscosity. 

Subprogram  FLIGHT  contains  the  differential  equations 
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of  flight  which  are  evaluated  by  sequential  calls  from 
RUNGE. 

Subprogram  CDRAG  develops  the  aerodynamic  coefficient 
of  drag  and  the  drag  increment  due  to  wings  (or  other 
aerodynamic  surfaces)  deployed  during  flight.  For 
accurate  simulation  of  a specific  system, tables  of  Mach 
number  and  coefficient  of  drag  are  read  by  MAIN  and 
transferred  thru  COMMON  to  CDRAG.  If  a tabular  drag 
coefficient  is  not  provided  as  input,  the  program  defaults 
to  an  internally  supplied  function. 

Subprogram  ACOEPS  develops  all  the  aerodynamic  co- 
efficients exclusive  of  drag  necessary  for  treating  spin- 
stabilized  projectiles.  An  internally-supplied  coeffi- 
cient is  used  whenever  tabular  data  for  that  coefficient 
is  not  entered  as  input.  The  coefficient  values  used 

in  default'  are  functional  fits  to  the  155  mm  (BAP)  pro- 
jectile T387. 

il-.  Options 

To  be  able  to  treat  a variety  of  applications  con- 
veniently a number  of  program  options  are  provided  — 
both  explicitly  as  option  switches  and  implicitly  as 
parameter  values  of  input  variables.  For  example,  a 
binary  switch  lOPTY  is  required  as  input.  If  this  para- 
meter value  is  omitted  or  set  to  0 (zero),  the  program 
bypasses  the  computation  of  normal-body  forces  required 
for  accurate  simulation  of  spin-stabilized  projectiles. 
Additionally,  the  program  does  not  expect  or  read  the 
input  parameters  needed  for  computing  projectile  spin 
and  yaw.  A considerable  savings  in  execution  time  is 
obtained  by  omitting  the  spin  option.  If  the  spin  op- 
tion is  desired,  lOPTY  is  sot  to  1 and  the  values  of 
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certain  (required)  variables  are  provided  to  Implement 
this  option. 

Another  major  option  Is  the  choice  of  purely  ballis- 
tic flight  or  flight  with  midcourse  glide.  This  option 
Is  exercised  by  setting  the  switch  IGLIDE  to  1,  The 
default  option  with  IGLIDE  set  to  0 (zero)  Is  ballistic 
fllght_.  If  IGLIDE  Is  set  to  1,  several  suboptions  are 
provided/required.  Since  midcourse  glide  or  body  alti- 
tude hold  applies  exclusively  (at  present)  to  the  Copper- 
head projectile,  default  tables  of  aerodynamic  coeffi- 
cients characterizing  this  system  In  glide  are  provided 
In  a BLOCK  DATA  subroutine.  To  override  these  entries  a 
switch  IDFUFO  Is  provided.  With  IDFUFO  set  to  1,  the 
program  expects  user-supplied  aero  data;  otherwise  Internal 
data  are  used.  Several  alternative  means  are  provided  for 
simulating  the  start  of  glide. 

In  suboption  1 a desired  glide  angle  (GLIDE)  Is  pro- 
vided and  the  program  calculates  the  body  angle  of  attack 
required  to  hold  body  attitude  after  reaching  the  desired 
glide  angle.  A value  of  the  tolerance,  e (EPSTHE)  rela- 
tive to  the  desired  glide  angle  Is  also  entered  for  this 
option.  See  Attachment  1 for  details.  Before  altitude 
hold  of  the  body  Is  Initiated  an  Initial  angle  of  attack 
Is  calculated  such  that  the  lift  force  equals  the  gravi- 
tational force  normal  to  the  velocity  vector  at  that  point 
In  the  flight.  This  option  Is  useful  when  timer  settings 
are  unavailable. 

In  suboption  2 a time  TENABL  (In  secs)  Is  provided 
and  GLIDE  Is  set  to  -90  degrees.  With  this  option  the 
calculations  for  angle  of  attack  to  hold  body  attitude 
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start  at  TENABL.  A value  of  TENABL  equal  to  that  at 
which  the  velocity  angle  0 = GLIDE  produces  a trajectory 
Identical  to  that  of  suboption  1,  Both  suboptions  1 and 
2 are  exercised  by  omitting  (or  setting  to  zero)  the 
value  of  MMCSW. 

Suboption  3 requires  the  program  user  to  provide  time 
called  TMOKMC  In  the  program.  This  Is  the  time  from 
launch  at  which  the  Internal  (MMC)  timer  sequence  starts. 
This  option  Is  exercised  by  setting  the  switch  MMCSW  to 
unity  and  setting  GLIDE  to  -90  degrees.  Although  Internal 
degrees  of  freedom  within  the  projectile  are  not  simulated 
with  a three-degree-of-freedom  model.  It  Is  Important  to 
properly  treat  their  kinematic  effect.  In  contrast  to 
other  suboptions,  suboption  3 faithfully  simulates  the 
effect  of  the  events  which  occur  after  T^.  These  events 
are  tabulated  below: 


SEQUENCE  OF  EVENTS  IN  COPPERHEAD  TIMER 


Event 

No. 

Time 

(s) 

Description 

1. 

T 

o 

timer  sequence  starts  (main  power  at  + 30v 

2. 

T 

o 

+ 

2 

roll  control  starts 

3. 

T 

o 

+ 

3 

seeker  gyro  splnup  Initiated 

4. 

T 

o 

+ 

attitude-hold  enable  switch  is 

closed 

5. 

T 

o 

+ 

5 

start  to  extend  wings,  apply  g 
free  gyro 

bias,  and 

As  far  as  projectile  kinematics  are  concerned,  the 
projectile  remains  ballistic  until'  event  number  4.  When 
event  number  occurs  at  T^  + ^(s),  a control  surface  de- 
flection In  pitch,  6p,  Is  calculated  continuously  and 
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applied  throughout  the  Interval  between  events  4 and  5. 

This  algorithm  bases  the  calculation  of  6 upon  the 
equilibrium  behavior  of  the  projectile  airframe  and  auto- 
between  events  4 and  5«  This  behavior  Is  approxi- 
mately one  In  which  pitch  deflection  of  the  fins  and  the 
associated  trim  angle  of  attack  are  proportional  to  the 
average  turning  rate  of  the  velocity  vector,  e . Thus, 

<5p  = K e (K,  a constant). 

The  turning  rate  6 Is  calculated  continuously  by 

• • ••  • •• 

e = (xy  - yx)/v  . 

An  average  value  of  e Is  calculated  by  exponential  smoothing. 
This,  of  course,  assumes  that  this  system  exhibits  nearly 
first-order  dynamics  between  the  events  4 and  5. 

Specifically, 

-i- 

1 + TS  ^ 


with  T a time  constant  (0.2  sec)  and  s the  Laplace  differ- 
ential operator.  Using  exponential  smoothing,  this  equa- 
tion Is  solved  Implicitly  by  the  recursive  procedure: 


with 


®1  ° * '^1®! 


0 = e 

o 


where  h Is  the  Integration  time  step,  l.e., 

h = ^1-1  + ^ • 

For  accurate  simulation  of  this  portion  of  the  system 
h £ 0 . 1 sec . 

Having  calculated  6^,  the  trim  angle  of  attack  (a^)ls  ob- 
tained by  linear  Interpolation  In  the  following  table. 
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TRIM  ATTACK*  VERSUS  MACH  NUMBER  AND  CONTROL  DEFLECTION 


Entries 

are 

(deg) 

Mach 

6 

(deg) 

Number 

0 

5 

10 

15 

0.5 

0 

7.4 

12.1 

17.7 

0.8 

0 

7-2 

11.8 

17-2 

0.9 

0 

7-4 

11.6 

17.1 

1.0 

0 

7-1 

11-3 

17.1 

With  all  program  options  the  normal  body  force  is 
calculated  from  the  normal  force  coefficient  at  trim. 

This  coefficient  is  a function  of  Mach  number,  M,  and  trim 
attack,  , 1 . e. , 

Values  'Of  Cj^  are  obtained  by  two-way  linear  interpolation 
in  a table  as  explained  in  Attachment  1. 


* Reference:  Wind  Tunnel  Data  Analysis  of  3/^  Scale  Model 
for  the  XM712  Frojectile,  27  Feb  76,  pp.  72-75- 
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5.  Input  Variables 


Program  data  input  is  supplied  by  pxmched  cards 
whose  content  is  described  below  in  sequential  order. 


CARD 

PORMAT 

COLS 

VARIABLE 

CONTENT 

1 

20A1; 

1-80 

TITLE 

description  of  aero- 
dynamic data,  option- 
ally supplied 

2 

'212 

1-2 

NTBL 

nxamber  of  entries  in 
drag  table 

3 - 

NARTBL 

ntunber  of  entries  in 
each  aero  coefficient 
table.  A zero  or 
blank  exercises  de- 
fault . 

2a 

8P10.0 

1 - 10 

11  - 20 

• 

• 

XMTBL(l) 

XMTBL(2) 

first  entry  in  table 
of  Mach  numbers  asso- 
ciated with  drag. 
Points  entered  low 
to  high. 

71  - 8o 

XMTBL(8) 

2b 

8P10.0 

1-10 

XMTBL(9) 

• 

continuation  of  Mach 
no.  table 

2c 


2d 


8P10.0 


8P10.0 


1-10 

71  - 8o 


XMTBL(NTBL)  last  entry 

CDTBL(l)  first  entry  in  coef. 
of  drag  table 

CDTBL(8) 


1-10 


CDTBL(9)  continuation  of 
drag  table 


CDTBL(NTBL)  last  entry  in 
drag  table 

2e  8P10.0  1-10  COTBL(l)  first  entry  in  table 

of  drag  increment  due 
to  wings 

2f  8P10.0  1-10  C0TBL(9)  continuation  of  table 

of  drag  increments 

COTBL(NTBL)  last  entry  in  drag 
increment  table 
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CARD 


FORMAT 


COLS 


VARIABLE  CONTENT 


If  NARTBL  ^ 0,  additional  cards  must  be  supplied  for  each  of 
the  aerodynamic  coefficients  listed  below. 

TMACH(I)  table  of  Mach  num- 
bers used  for  the 
additional  aero  data. 
Each  table  must  have 
NARTBL  entries  with 
blank  cards  supplied 
for  missing  data. 

All  aerodynamic 
coefficient  values 
In  the  following 
tables  are  the 
balllstlcians ’ values 
based  upon  caliber 
squared . 

TKA(I)  table  of  spin  damp- 

ing moment 
coefficient 


311, 7a,  1 

3F10.0 


TKDYAW(I) 

table  of  yaw  drag 
coefficient 

TKL(I) 

table  of  lift  force 
(derivative)  coeffi- 
cient 

TKM(I) 

table  of  overturning 
moment  (derivative) 
coefficient 

TCP(I) 

table  of  center  of 
pressure  (In  cali- 
bers aft  of  nose) 

TKF ( I ) 

table  of  Magnus  force 
coefficient 

TKT(I) 

table  of  Magnus  mo- 
ment coefficient 

TKH(I) 

table  of  damping  mo- 
ment coefficient 

TKS(I) 

' table  of  pitching 
force  coefficient 

IGLIDE 

switch  for  glide 
option 
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CARD 

FORMAT 

COLS 

VARIABLE 

CONTENT 

2 

IDFUFO 

switch  to  Input 
aero  data  for  glide 

3 

MMCSW 

switch  for  glide 
suboption  3 

11  - 20 

TMOMMC 

timer  setting,  (sec), 
for  suboption  3 

21  - 30 

GLIDE 

desired  glide  angle, 
(deg).  Entered  as  a 
negative  value. 

31  - 40 

EPSTHE 

tolerance  In  the  de- 
sired glide  angle, 
(deg) 

If  IDFUFO  = 1,  additional  cards  must  be  supplied  for  each  of 
the  aerodynamic  coefficients  listed  below. 

3a 

4F10.0 

TFMACH(I) 

table  of  Mach  numbers 
used  to  Interpolate 
for  normal  body  force 
during  the  midcourse 
glide  phase  of  flight 

3b 

4F10.4 

ALTRMM(I) 

table  of  values  of 
maximum  trim  angle 
of  attack,  (deg) 

3c 

4F10.4 

AALTRM(I) 

argument  values  of 
angle  of  attack, 
(deg).  In  table  of 
normal  force  coeffi- 
cient versus  Mach 
number  and  trim 
angle  of  attack 

3d 

4F10. 4 

TCN(I,J) 

table  of  normal 
force  coefficient 
with  I Indexed  over 
Mach  numbers  and  J 
Indexed  over  trim 
angles  of  attack, 
(deg) 

4 

20A4 

1-80 

TITLE 

description  of  run 
set 

5 

8F10.0 

1-10 

D 

caliber  or  refer- 
ence diameter,  (mm) 

11  - 20 

EMO 

Initial  projectile 

mass,  (lb) 
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CARD  FORMAT  COLS 

21  - 30 

31  - ^0 
Hi  - 50 

51  - 60 

6l  - 70 

71  - 80 

6 7P10.0,  1-10 

313 

11  - 20 
21  - 30 

31  - ^0 

Hi  - 50 
51  - 60 

61  - 70 

71  - 73 
JH  - 76 
77  - 80 


VARIABLE 

CONTENT 

EMB 

burnt  or  final  mass, 
(lb) 

PC 

nominal  thrust,  (lb) 

SPI 

specific  Impulse, 

( sec ) 

DELT 

rise  time  of  thrust 
to  nominal  value, 

( sec ) 

VO 

Initial  velocity, 
(f/s) 

VWP 

velocity  of  constant 
headwind,  (f/s) 

HO 

altitude  ASL  of 
launch  point,  (f) 

HTERM 

altitude  ASL  of 
Impact  point,  (f) 

PPCTR 

form  factor  or  ratio 
of  ref.  area  used  In 
est.  aero,  coefs.  to 
ref.  area  used  In 
Input  data 

QEO 

Initial  value  of 
quadrant  elevation 

In  parameter  set, 
(deg) 

DQE 

Increment  In  QE, 

(deg) 

STEP 

time  step  supplied 
for  numerical  Inte- 
gration, (sec) 

TM 

time  programmed  for 
start  of  thrust, 

( sec ) 

NQE 

number  of  steps  In 

QE  desired 

NPRINT 

number  of  time  steps 
per  printout 

lOPTY 

switch  for  spin 
option 
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CARD  FORMAT  COLS 

7 4F10.0  1-10 

11  - 20 

21  - 30 

31-^0 

If  lOPTY  Is  zero,  following  cards 

7a  6F10.0  1-10 

- 12 

11  - 20 

21  - 30 
31  - ^0 
4l  - 50 

51  - 60 
61  - 62 

7b  6F10.0  1-10 

11  - 20 


VARIABLE 

CONTENT 

CADENS 

correction  factor 
for  air  density 
relative  to  NASA 
std.  atmos. 

VCW 

velocity  of  cross- 
wind  from  right  to 
left  facing  down- 
range,  (f/s) 

TENABL 

time  at  which  mid- 
course glide  Is  en- 
abled, (sec) 

THID 

thrust-induced-drag 
factor  applied  dur- 
ing burning.  Set 
to  unity  for  stan- 
dard case  . 

are  not 

required. 

SPINO 

Initial  spin, 

(rad/sec ) 

XCG 

position  of  the  cen- 
ter of  gravity  aft 
of  nose,  (cal) 

CLONG 

length  of  projectile, 
(cal) 

AMOM 

axial  moment  of  Iner- 
tia, (kg/m2) 

BMOM 

trasverse  moment  of 
Inertia  thru  eg, 
(kg/m2) 

WTAREA 

ratio  of  wetted  area 
to  reference  area 

ISEP 

switch  for  separate 
computation  of  skin 
friction  drag 

VHXF 

velocity  of  the 
launcher  In  the  X- 
dlrectlon  In  Iner- 
tial space,  (f/s) 

VHYF 

velocity  of  the 
launcher  In  the  Y- 
(vertlcal)  direction 
In  Inertial  space, 
(f/s) 
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CARD 


FORMAT 


COLS 

VARIABLE 

21  - 

30 

RCW 

1 

H 

ko 

XCW 

CONTENT 

scale  coefficient  in 
rangewise  windshear 

scale  coefficient  in 
crosstrack  windshear 
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6.  Output  Variables 

The  definition  of  program  outputs  and  echoed  inputs  is 
listed  below. 


FORTRAN 
NAME 


UNITS 


DESCRIPTION 


TENABL 


sec 


time  at  which  midcourse  glide 
is  enabled 


THID 
C A DENS 
PPCTR 


nondimens ional  thrust  induced  drag  factor 

(nominally  unity) 

nondimens ional  correction  factor  to  air  den- 
sity (nominally  unity) 

nondimens ional  factor  used  to  adjust  drag 

coefficient  for  non-standard 
conditions,  for  example,  when 
aero,  coefs.  were  developed 
for  a different  reference  area 
than  that  used  for  a run 


VO 

f/s 

EMO 

lb 

m 

EMB 

lb 

m 

D 

mm 

QE 

deg 

DT 

sec 

PC 

SPI 

Ib-sec/lb 

VWP 

f/s 

HO 

m 

HTERM 

m 

VCW 

f/s 

initial  (muzzle)  velocity 

initial  mass  of  projectile 

final  or  burnt  mass  of  pro- 
jectile 

caliber  or  reference  diameter 
of  projectile 

quadrant  elevation  or  launch 

with  respect  to  launcher 

integration  time  step 

nominal,  constant  level  of 
thrust 

specific  impulse  of  rocket 

velocity  of  headwind 

initial  altitude  above  mean 
sea  level  (KSL) 

terminal  or  target  altitude 
above  MSL 

velocity  of  crosswind  (from 
right  to  left  facing  downrange) 
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FORTRAN 

NAME 

UNITS 

DESCRIPTION 

RWC 

nondimensional 

scale  factor  for  windshear 
function  for  headwinds 

XWC 

nondimens i onal 

scale  factor  for  windshear 
function  for  crosswinds 

T 

sec 

time  after  lavinch 

X 

m 

rangewise  pro;jectile  coordinate 

HAG 

m 

height  above  ground  Impact 

Z 

m 

crosstrack  or  deflection  co- 
ordinate 

XD 

m/s 

X-component  of  pro;jectilo 
velocity  (Earth  coordinates) 

YD 

m/s 

Y-component  of  pro;jectile 
velocity 

V 

m/s 

pro;jectile  speed  in  the  X-Y-Z 

inertial  frame 

CMACH 

nondimensional 

Mach  nvimber  of  projectile 

RESTS 

lb 

drag  force 

THETA 

M 

deg 

attitude  of  velocity  vector  of 
projectile  in  vertical  plane 

YAWDEG 

deg 

angle  of  attack  (generalized 
yaw)  of  the  projectile  Zero 

is  returned  for  lOPTY  = 0. 

SPIN 

rad/s 

projectile  spin  Zero  is 

returned  for  lOPTY  = 0. 

DMY 

Dummy  variable  is  provided  for 
convenience  of  user. 

SUPVEL 

m/s 

maximum  speed  of  the  projectile 

RANGE 

m and 

range  of  projectile  relative  to 

nautical  miles 

t surface  of  Earth 

SUPALT 

m 

maximum  altitude  relative  to  MSL 

PITCH 

deg 

pitch-plane  angle  of  attack 

BODYA 

deg 

pitch-plane  body  attitude  in 
Earth  coordinates 
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FLOW  CHART  FOR  EXTERIOR 


BALLISTICS  OF  BOOSTED  ROCKETS 


1.  Read  in  constants  and  params . 

2.  Compute  auxiliary  constants  and  redimension  input  data. 

3.  Set  constants  for  QE  - loop. 

4.  Start  QE  - loop. 

5.  Print  and  label  variables. 

6.  Set  initial  conditions  for  time  loop,  e.g.,  ISW  = 0 & 
t,x,y,x,y . 

7.  Call  RUNGE  1. 

8.  Initialize  NPL0T. 

9.  Call  subroutine  of  flight  equations,  FLIGHT,  for 
evaluation  of  initial  conditions. 

10.  Initialize  a counter  LINE  for  counting  lines  on  a page 
for  proper  labeling  of  output  at  top  of  each  page. 

11.  Skip  to  16  and  print  initial  conditions. 

12.  Start  time  loop. 

13.  Start  print  loop. 

14.  Test  if  burning  has  begun;  if  so  set  T0  = TIME  and 
bypass  14  in  future. 

15.  Call  RUNGE  2 for  solving  flight  equations. 

16.  Count  lines  printed  and  check  if  time  to  skip  to  new 
page  and  label. 

17.  Print  time  and  dependent  variables. 
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18.  Inci’einent  NPL0T  and  store  position  of  projectile  in 
XPL0T  and  YPL0T . Save  projectile  position. 

19.  Check  if  time  to  stop  integrating.  Stop  when  y goes 
thru  YTERM.  If  y>YTERM,  return  to  step  12. 

20.  If  y^YTERM,  use  past  values  of  dependent  variables  to 
interpolate  linearly  for  their  values  at  YTERM. 

21.  Print  out  final  values  of  dependent  variables. 

22.  Call  the  plotting  subroutine  PPL0T  and  ask  for  a plot 
of  the  trajectory,  represented  by  the  x,y  pairs  saved 
in  XPL0T  and  YPL0T . Label  the  plot  with  the  run 
description,  the  quadrant  elevation,  initial  velocity 
and  nominal  thrust  level. 

23.  Return  to  1 for  another  parameter  set,  otherwise 

24.  S top . 
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Flow  Chart  for  Subroutine  FLIGHT 


1.  Teat  if  ISV;  = 1. 

2.  a.  If  ISW  = 0,  omit  thrust  terms. 

b.  If  ISW  =1,  include  thrust  terms. 

3.  a.  If  ISW  = 0,  use  initial  mass 

EM  = EMO. 

b.  If  ISW  = 1,  call  subroutine  BURN.  CALL  BURN(TIME, 
XMASS,  THRUST). 

This  generates  the  projectile  mass  in  povinda  mass  and  thrust 
in  pounds  force.  Conversion' of  units  occurs  in  the  differen- 
tial equations. 

k.  CALL  SOUND  (Y,  A,  G,  RHO ) . 

This  generates  speed  of  sound  (A),  gravity  (G),  and  air  den- 
sity (RHO) . 

5*  Relative  wind  speed  computed:  VREL  = SQRT( (XD0T+VW)*v2 
+ YD0T-"--::-2  + ( ZDOT+VCW )*-::- 2) 

6.  CMACH  = VREL/A 

7.  a.  If  lOPTY  = 0,  omit  spin  and  normal  force  calculations. 

b.  If  lOPTY  = 1,  call  ACOEFS  and  calculate  spin  and 
normal  force  derivatives. 

8.  a.  If  ISEP  = 0,  set  FRICT  to  zero. 

b.  If  ISEP  = 1,  compute  akin  friction  drag  (FRICT) . 

The  subroutine  CDRAG(CMACH,  DRAG,  CAD)  generates  the  un- 
corrected coefficient  of  drag,  DRAG,  and  the  drag  increment 
due  to  wings,  CAD. 

9.  COFDRG  = FFCTR  DRAG  + XKDYAW  YAWSQ  + FRICT  + CAD. 

10.  a.  If  IGLIDE  = 1 and  NABLE  = 1,  compute  normal  force 

terms  in  GLIDE  mode. 

b.  Otherwise,  omit  GLIDE  calculations. 

11.  Solve  differential  equations  for  second  derivatives  of 
state  variables. 

12.  Return. 
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ATTACHMENT  1 


Algorithm  for  Attitude-Hold  Logic  In  Copperhead  3D0F  Simulation 

For  simulations  in  which  a preset-time  option  is  not  used,  the 
program  determines  the  time  for  commencement  of  attitude  hold.  Initially 
the  glide  angle  0 will  approach  the  desired  glide  angle  0^  algebraically 
from  above.  (0o  is  negative)  When  0 - 0o  1 £ (e  positive),  an  initial 
angle  of  attack,  a^,  is  computed  such  that  the  associated  lift  component 
equals  the  component  of  gravity  normal  to  the  velocity  vector,  ie,  such 
that 

F,  “ M g cos  0 
L p 

with 


and 


F_  ” F„  cos  a 
L N o 


with  projectile  mass  M and  reference  area  A and  dynamic  pressure  q. 

Since  is  small—typically  less  than  10® --an  Iterative  procedure 

is  employed  in  which  the  first  iteration  assumes  that  F^^  * F^^  so  that 

(1)  M g cos  0 
— 

The  value  of  is  obtained  by  interpolation  in  the  tabular  function 

o 

C„(M,  a ) with  a the  trim  angle  of  attack. 

N t t 


Thus, 


(1)  * (1) 
C„(0)  = C^,(M*,  ^), 


'N'“'  N'  ' o 

where  M*  is  the  local  Mach  number.  The  interpolation  procedure  is  described 
below. 

Then,  form  the  second  iterate  for  the  normal  force  coefficient: 

(2)  M g cos  0 
C„(0)  = 


N 


Aq  cos  a 


(1) 
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The  value  of  obtained  by  requiring  that  taken 

as  the  initial  trim  angle  of  attack  at  the  start  of  attitude  hold. 

Interpolation  Procedure 

To  calculate  a^,  interpolate  on  Mach  number  in  the  ct^)  table 

obtaining  a^)  at  the  local  Mach  number  M*  for  values  of 

Ct^  “ {O,  5,  10,  15  (deg)}.  Then,  is  obtained  by  linearly  interpolating 
with  Cjj(O)  argument  in  this  table. 

The  value  of  initial  body  attitude  is  given  by 

6^(0)  - e + o„. 

For  subsequent  calculations  during  the  attitude-hold  trajectory, 

the  value  of  angle  of  attack,  a,  is  calculated  which  preserves  body 

attitude,  ie,  for  which  9,  * 9,  (0) . 

D b 

Thus, 

a * 0.(0)  - 9. 

D 

If  the  above  value  of  a satisfies 

ot  < a^^(M) , the  maximum  trim  angle  at  M the  instantaneous  value 
of  Mach  number,  the  lift  force  is  computed  as 

F * C (a,  M)  Aq  cos  a 

Lt  N 

with  obtained  by  two-way  interpolation*  Otherwise,  a is  limited  to 

a (M)  and  lift  is  calculated  as 
tm 

’’l  ' “tm- 

Induced  drag  due  to  lift  is  calculated  as 
Fjjj  = - M)  Aq  sin  a . 
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TABLE  1.  MAXIMUM  TRIM  ANGLE  OF  ATTACK 
VERSUS  MACH  NUMBER  (AT  6 - 12®) 


Mach 

Number 

(max) 

(deg) 

0.5 

14.3 

0.8 

14.0 

0.9 

13.8 

1.0 

13.6 

TABLE  2.  NORMAL  FORCE  COEFFICIENT 
VERSUS  MACH  NUMBER  AT  TRIM  ATTACK 


Mach 

Number 

(deg) 

0 

5 

10 

15 

0.5 

0 

1.1 

2.1 

3.1 

0.8 

0 

1.2 

2.1 

2.9 

0.9 

0 

1.3 

2.2 

3.0 

1.0 

0 

1.3 

2.5 

3.7 

Next  page  is  blank. 
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ANNEX  2 


SOURCE  PROGRAM  LISTING  WITH  EXAMPLE 
FOR  THE  COMPUTER  PROGRAM: 
EXTERIOR  BALLISTICS  OF 
BOOSTED  ROCKETS  (EXBAL) 
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no  o o n o n 


G LEVEL  21 


MAIN 


DATE  = 77143 


17/22/29 


EXTEf^IOR  ballistics  OF  BOOSTED  ROCKETS 
A THPEE-OEGRFE-OF-FREEDOM  MODEL  APPLICABLE  WHERE 
trailing  or  FOLLOWir>JG  FltHAVIOR  CAN  Bl  ASSUMED 


real  RS (401 ) »TS (401 ) » vs (401 ) ,CJRT ( 1 1 ) ,CJVT (11) 

OIMEMSION  TITLE(20) ,(I(I?) ,wP(48) ,XmI-L(]?) »CDTBL(I2) *C0TBL(I2) 

1 » Tl'iACH  (11),  Ti:A  (ID,  TKDYAW  (11),  TKL  ( 1 1 > ,TKM  (II)  ,TKF  (11), 

2 TKTDI)  ,TKH(n)  ,TKS(11)  ,TCP(II)  ,TFMACH(4)  ,ALTRMM(4)  ,AaLTRM(4)  , 

3 TCN(4,4) ,WTCN(4) 

DIMENSION  TC4D1  (4)  ,TCAL)2(4)  ,TDEL(4)  ,TaT(4,4)  ,WTAT(4) 

INTEGER  *2  Ch AR  ( 1 )/•<>•  / 

DATA  RF/6.37BE6/  ,OMF.GA/0.72915E-4/ 
data  DTORAD/S7. 29578/ 

DATA  lMCONG/0.2/ 

ASSIGN  CONSTANTS  NEEDED  BY  DIFFERENTIAL  EQUATIONS  TO  COMMON 

COMMON  EM0,E:-U'},SPI  ,FC,HRATE,nELT,TO»Te,  I SW  , V , THE  T A , FFCTR  , CaLSQ  , 

1 VW,vCW,ALT,R,  IENn,CNACH,REY''JLD,RESlS,CAL,OLONG,  I OPT  Y , YAW  , aMOM  , 

2 i-;MOf,PSl  ,wTAREA,  ISEP,MAr-LF.,  IGLIDE,  DlNG,TENAbL,  AVTH0,CDEFl, 

3 M'^CSW,  TmOMMC 
CO^M()N  /SRCOi./lM 

CO^'MO.i/CuFCOM/XCG,SMARG,EM,  THID,PRnU  , aL  TR  I M , CN  A TRM  , GL  I DE  , EPSTHE 

1 ,Sr  aF  AC,  YA-vNU,  TKA,  TKUYAW  ,TKL»  TKm,  TKF  , TKT  , TKH  , TK  S , TCP  , 1 MACH , 

2 NaRTKL  , TFmACH,AI.TRMM,  AALTKM,TCN,vD!rN.TCADI  , TC  AD2  , TUEL  , TAT  , WT  AT 
C OrMON/ w I NCOM/R  w C , X w C 

COMM(iN/l)RGC0'i/  XMTyL,CDTBL,COTBL,NTHL 

COMMON/ SNOCOr'./CADENS 
C 

c*<nnt  Tables  of  aerodynamic  coefficienis  is  param.  set  input  set  -i. 
c<nnnt  IF  Parameters  ntpl  and  nartbl  are  poih  zero,  endogenous 

C«n><nt  FDI'.'CTIONAL  fits  TO  TtlE  AERODYNAMIC  TABLES  (WITH  ThE  T3H7  FORM) 

c<nnn»  ,^ii_L  be  usf‘0.  see  subroutine  acolfs. 

Cinn>*  IF  OfsLY  NARIBL  is  ZERO,  ThE  ZF.RO-LIFT  DRAG  TABLE  IS  REQUIRED 
C**o«  WITHOUT  REQUIRING  TAbLES  FOR  THE  OTHER  AERO  COEFFICIENTS. 

C«<»**  IF  CERTAIN  AERO  COEFFICIENTS  ARE  UEEir,'En  (KNOWN),  THESE 
C<nt«*  CAN  bE  read  with  THE  OTHERS  LEFT  bLA  K.  THE  PROGRAM  wiLL 
C<nnn>  USE  the  TaBULATEN  COf.FF  I C I ENTS  aND  DEFAULT  TO  THE  ENDOGENOUS 
F'UNCTIONS  FOR  THOSE  ENTERED  AS  ZERO* 

c o=projectile  caliber,  millimeters 

C EMO= initial  projectile  mass,  lbm 

c EmH=pURNI  mass,  lhm 

C FC=T.'OMINaL  Thrust  level,  Lf;F 

c spi=specific  impulse  of  rocket  propellant, 

C F.RaTF.=PROPELLANT  riURNING  RATE,  LBM/BEc 

C DELT=ThRUST  RISE  TIME,  SEC  ’ . ' 

C T0=IGNITI0N  TIME  FOR  ROCKET  MOTOR  ENDOGENOUS 

c IN  Subroutine  *hurn'  the  thrust  decay  time  is  assumed 

c equal  to  the  Thrust  rise  time,  a typical  value  = o.i  sec. 

C T8=eFFECTIVE  ‘'URNliNG  INTERVAL*  SEC  ENDOGENOUS 

C ISw=  A SWITCH  signaling  COMMENCEMENT  OF  BURNING  ENOO. 

C IEND=A  SWITCH  SIGNALING  ENU  OF'  BURNING  ENDO. 


parameter 


LBF/LBM/SEC 

ENDOGENOUS 


INPUT 

INPUT 

INPUT 

INPUT 

INPUT 


OOOOOlOO 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

ooooiooo 
00001 100 
00001200 
00001300 
00001400 
00001500 
00001600 
00001700 
OOOOIBOO 
00001900 
00002000 
00002100 
00002200 
00002300 
00002400 
00002500 
00002600 
00002700 
00002800 
00002900 
00003000 
00003100 
00003200 
00003300 
00003400 
00003500 
00003600 
00003700 
0000-800 
00003900 
00004000 
00004100 
00004200 
00004300 
00004400 
00004500 


VaRI ABLE00004600 
INPUT  6 00004700 
VARI ABLE00004800 
00004900 
00005000 
VARIABLE00005100 
VARIABLE00005200 
VARIABLE00005300 
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ooo  onoooo  onoonnnooooooooooononoonnnono  onono 


G LEVEL 
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MAIN 


DATE  = 77143 


; V=PKOJECTILE  VELOCITY,  M/SEC  ENOO, 

VO=MliZZLE  VELOCITY  OF  TriL  PROJECTILE,  FT/SEC. 

theta=attituoe  of  prujectile,  deg 

CALSi:=CALIt->EW  SOUArVEl),  N*<>2 

VW=VfLOCITY  or  HFADNIND,  M/SEC  (REaU  IN  IN  FT/SEC) 

Vwr  = VELOCITY  <JF  HEAD'^INO  IfJ  FT/SEC 
IEY033I  CO'-U'ENTS  DELETED 

continue 

alt=true  altituof  afiove  sfa  level,  m 

R=RAOIUS  FROM  CENTER  OF  EARTH  TO  PROJECTILE  , M 
RE=N0MINAL  RADIUS  OF  TllE  EARTH  AT  THE  EQUATOR,  M 
OMEGa=ANGUL AR  VELOCITY  OF  THE  EARTH,  PAD/SEC 
IOPTy-  a S.VITCH  INDICATING  CHOICE  OE  YAW  OPTION. 

I0PTY=  I PRODUClS  COHPUIaTiON  of  yaw  of  repose  FOR  spinning  PROJEC00007600 
I0PTY=  0 SIGNIFIES  A TRAILING  PROJECThE  WITHOUT  SPIN,  FOR  00007700 


ENOO. 

ENDO. 

ENDO. 


ENDO. 

ENDO. 


VARIABLE00''05400 
INPUT  7 00005500 
VARI  A8LE00005600 
VARIABLE00005700 
VARI ABLE00005800 
INPUT  8 00005900 
« * « * « * « « , 

00007000 
VARIABLE00007I00 
VARI ABLE00007E00 
CONSTANT00007300 
CONSTANT00007400 
INPUT  1800007500 


This  option  the  following  in^^uts  are  unnecessary. 

SPIN(,i=INITIaL  spin,  rad/sec 

XCG=P()SI  T lUN  OF  center  OF  GRAVITY  aFT  OF  NOSE,  CALIBERS 
XCP=F>OSITION  of  center  of  pressure  aft  of  nose,  cal  ENDO. 
CLONO=PROJECT ILE  LENGHT  IN  CALIBERS* 

AM0M=l0NGI TUDIUAL  MOMENT  OF  INERIIA  OF  THE  PROJECTILE,  KG«M«*2 


00007800 
INPUT  1900007900 
INPUT  2000008000 
VARI ABlEOOOOBI 00 
INPUT  2100008200 
00008300 


INPUT  2200008400 

BMUF'rTRANSVtPSF  MOM.ENT  OF  INERTIA  OF  THE  PROJECTILE,  KG*MHn*2  00008500 

INPUT  2300008600 

VCw=VELOCITY  OF  CROSSwiND  FROM  RIGHT  LOOKING  DOWNRaNGE (READ  IN  FT/00008700 

INPUT  2400008800 

wtaRe A=wETTEr>  apfa  ratio  used  in  computing  skin 

FRICTION  I^PaG 

ISEP=A  Switch  indicating  choice  of  separate 
computation  of  skin  FRICIION  drag. 


00008900 
INPUT  2500009000 
00009100 
INPUT  2600809200 


= 1 IF  FRICTION  DRAG  IS  CU'^PJTLD  SEBARATELY  AITD  ADDED  TO  FORM  DRAG00009300 


= ^ IF  friction  drag  is  included  in  urag  function. 
SMAR(R=PROJlCTILF  static  MAhGIN,  CAL* 

PSI=rNGULAR  URIENTATIO^J  t)F  YAW  VECTOk 
SRNG=SL.ANr  RANGE  TO  PROJLCTIL'^  POSITION,  M 


VHXF=VELUCTTY  Of  LAUNCHER 
VHYF =VELOCITY  OE  LAUNCHER 
CONI INUE 

RWC  IS  HEADWIND  COEF. 

XWC  IS  CROSSWIND  COEF. 

PSI  may  be  computed  by  removing  'C 
IN  SURROUTIME  flight. 


IN  RANGEwISE  DIRECTION 
IN  vertical  direction. 


.00009400 
ENDO.  variable;oooo95oo 
ENDO.  VARIABLE00009600 
ENDO.  VARI ABLE00009700 
(FT/SEC)  00009800 

(FT/SEC)  00009900 

OOOIOOOO 
OOOlOlOO 
00010200 

S FROM  comment  CARDS  00010300 

00010400 

CAOEilS  IS  THE  CORRECTION  FACTOR  FOR  AIR  DENSITY  RELATIVE  TO  STANDAOOOI 0500 

00010600 

external  FLIGHT 

equivalence  (U(1),X),(U(2),Y),(U(3),Z),(U(4),SPIN), 

1 (U(5) ,XD)  , (U(6)  ,YD)  , (U(7)  ,ZD)  , (U(8)  ,SPIND)  , 

2 (U(9)  ,XDU) , (U(in) ,YDO)  , (U(1I ) .ZOD)  » (U(I2)  ,SDn) 

READ  In  run  Description,  constants  In  flight  equations  and 
initial  Conditions. 


read  (5,256,END=30)  T I TLE , NTBL , NARTBL 
256  format (20A4/2I2) 


00010700 
OOOIOfiOO 
00010900 
0001 1000 
OOOIIIOO 
000II200 
0001 1300 
0001  1400 
0001 1500 
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IF (MTBL.EO.O)  GO  TO  255 
read  (5»250)  (XMTHL ( 1 ) * 1 = 1 *NTBL) 
read  (5»250)  (CDTBL ( I ) » 1=1 ♦NTBL) 
read  (5»250)  (COTBL ( I ) ♦ 1 = 1 ♦NTHL) 

250  format (bFlO.O) 

WRITE  (6,252)  title 

252  FOR^'AT  ( 1 h 1 ,20  AA/ 1 HO  , 1 OH  MACH  NO,lOH 

DO  253  1 = 1 ,NTBL 
WRITE  (6,254)  XMTilL  ( I ) ,CDTbL  ( I ) ,COTtiL  ( 1 ) 

253  COfiTiNUt; 

IF  (NaRTBL.EO.O)  go  to  255 


COEF  DRAG,10H  drag  INCR) 


read 

(5f250) 

(Tmach ( I ) , 1=1 ,nartbl) 

read 

(5f25U) 

(TKA  ( I)  , 1 = 1 ,NaRT8L) 

read 

(5f25n) 

(TKC'YAW  ( 1 ) , 1 = 1 ,NARThL) 

read 

(5f 250) 

(TKL  ( I ) , 1 = 1 ,NaRTBL) 

read 

(5^250) 

(TKm ( I ) , 1 = 1 ,NaRTBL) 

READ 

(5f250) 

(TCP ( I ) ,1  = 1 ,NARTBL) 

read 

(5f 250) 

(TKF  ( I ) , 1 = 1 ,NaRTBL) 

read 

(5t250) 

( TKT ( I ) , 1 = 1 ,NaRTBL) 

read 

(Sf250) 

(TKr)  ( I ) , 1 = 1 ,NaRTBL) 

read 

(5f250) 

( TKS ( I ) , 1 = 1 ,NaRTHL) 

254  FOR'^aT  ( 1H  ,3F  10.4) 

WRI  TF 

(6,272) 

N0,8X,2HKa,5x,5hKDYAW, 
CP,  CAL) 


272  F0R»’A  T ( iHO  , 1 OH  MACH 
1 HX  , 2)HKL  , HX  , 2hK f ‘ , 1 OH 
DO  25  7 I = l,r.'ARTHL 

T--’ACH(I)  ,TKA(I)  ,TKDYAW(i;  , TKL  ( I ) , TKM  ( I ) ,TCP(i) 
FOfxf'AidH  ferio-S) 

co^itinue: 

CON’TlNUt 


251 

257 

255 

C 

C 

C 

C <0- <}►•{>  •{> 

C * 

C o o <> 
C o 

101 

264 

1264 

1 


1265 


GLIDE  - FUFO  - GLIDE  - FUFO  - 
ATTITUDE  HOLD  LOGIC  FOR  NEARLY 


the  s.vitch  ioliof;  must 
the  S'.vnCH  lOFUFO  MUST 
BE  Pk'ovioed  For  glide 
THE  PARAMETER  '-"-'CSvm  IS 
OR  GLIDE  trajectory  IS 


GLIDE  - FUFO  - GLIDE  - FUFO 
CUDSTanT  GLIDE  ANGLE  APT? 


GLIDE 


00011600 

00011700 

00011800 

00011900 

00012000 

00012100 

00012200 

00012300 

00012400 

00012500 

00012600 

00012700 

00012800 

00012900 

00013000 

00013100 

00013200 

00013300 

00013400 

00013500 

00013600 

00013700 

00013800 

00013900 

00014000 

00014100 

00014200 

00014300 

00014400 

00014500 

00014600 

00014700 

00014800 

00014900 


BE  SET  TO  1 L .R  an  ATT  I TUOE-HOLD  TRaJECTORYOOO  1 5000 
RE  SET  TO  1 If  NORMAL  FORCE  AERO  IS  TO  00015100 

OTHERwISl,  uEFAULT  AERO  DATA  ARE  USED.  00015200 
A SaITCH  To  dE  SET  TO  1 WHEN  AN  A T T 1 TuDE-HOOOO 1 5300 
TO  f(E  SIMULaIEd  using  TIMER  DaTa  SUPPLIED  B00015400 


TnE  ‘‘ARTIn  marietta  CORp.  (MMC)  . 1,.  USING  THIS  OPTION 

The  time  PAwamfIFR  TfiOMMC  hEPRESENTG  the  time  in  secs  from 
LAU^.^CH  at  ,VHICf(  THE  30  VOLT  POwER  BECOMES  AVAILABLE. 

Suf-SFQUpTT  EVtNTS  IN  SEQUENCE  ARE  TRlATED  IN  SUBROUTINE  FLIGHT 

CU<\'T  1 ,vjUL 

REaU  ( 5,264,  End  = 30  ) IGLIl^E,  I DFUFO  , MMCSW  , TMOMMC  , GL  1 OE  , EPSTHE 
format (311 ,7x,1F  10.0) 

write  (6,1264)  IGL IDE , IDFUF 0 ,MMC5w , TMOMMC ,GL1 DE , EPSThE 
format  ( lHO,4X,6HinLIDE,4X,6HIDFUFO,5X,5HMMCSW  ,10H  TOMMC  (S) 
lOH  glide  (U)  ,11h  EPSTHE  ( D ) /3 ( 9X , 1 1 1 ) , 3F 1 0 . 4) 

IE  ( lUFUFO.NE. 1 ) GO  TO  1300 
read  (5,1265)  (TEMACH(I) ,1=1 ,4) 

Format (afio.o) 

read  (5,1265)  (ALTRMM ( 1 ) , 1=1 ,4) 


00015500 

00015600 

00015700 

00015800 

00015900 

00016000 

00016100 

00016200 

00016300 

00016400 

00016500 

00016600 

00016700 

00016800 
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READ  {5»1265)  { A ALTRM { I ) , I = 1 , 4 ) 

read  (5*1266)  ( (TCN { I * J)  * J=1 *4) , 1 = 1 *4) 

1266  format (4510.0) 

1300  continue 

IF  ( IGLIOE.UE. I)  GO  TO  1 
WRITE  (6,26ti) 

268  FORMAT  ( I)*0.8X. 'MACH  tviO  MAX  TRIM  (Q)') 

WRTTf  (6,269)  ( T FMACM  ( I ) * AL  TRm'1  ( I ) , I = I . 4 ) 

269  format (IN  ,2F 15.5) 
wRiTf.  (6,270) 

270  format  ( 1H0,10X, '^^ORMAL  FORCE  COEFS  AT  TRIM  ATTACK’/ 

1 IH0.8X,25HMACH  no  TRIM  ATTACK  (U):  * 1 (TO  , 14X  , lH5 , 1 3X  , 2H 1 0 , 1 3X  , 

2 2HI5) 

WRITE  (6,271)  (TFMACM(I) , (TCN(I,J) ,0=1,4) ,1=1,4) 

271  FORMaT(1)H  ,F15.5,4X,4F15.5) 

1 READ  (5,2,ENO  = 30)  T I TLE , U , E mo  , EM8  , FC , SP I , OELT , VO , VWF , HO  * HTERM , 

1 FFCTR,OF.O, DOE, STEP,  T M , NOE  , NPR  I NT  , lOPTY 

2 FORMAT (20A4/ttF10.0/FFlO.O,3I3) 

C<nnn*  Sv-ITCH  NATiLL  IS  SET  FROM  0 TO  1 AT  TIME  TENaBL 
256  PEaD  (5,260)  C AOENS , VC W , TEN ARE , TH I D 
260  forma  I (4F  10.(1) 

WkiTF  (6,262)  TClJAbL,  TfiIU,CADENS 
262  FORMAT  ( IHO,  l4fTENAMLF  TIMf.  = ,Fl0.4,5x, 

1 21(aHRUSr  IJOAO  FACTOR  = ,F10.4,5X,ic3HAIR  DENS  FACTOR  = ,F10.4) 
IF  ( lOPTY.rjF  . I ) 00  lO  25 

READ  26,SPINO,XCG,CLONG,AMOM,0mOm,wT  -.REA,  ISEP 
READ  26, VHXF,VHYF,RwC,XwC 

26  format (6F I 0.0, 12) 

GO  TO  27 

25  SPINO=0.0 
XCG=0. 

CLONG=0. 

AMOM=0. 

dMOM=0. 

WTARFA=0. 

ISEP=0 

vhxf=o.o 

VHYF  =0. 0 
RWC=0.0 
XrtC  = 0.() 

SmaR(.i.=  0 . 0 
staf  aC=0.0 
YAWNII=0. 0 
PRNU=0.0 
PSI=0.0 
OMY=0.0 

C START  OUAD-ELEV  LOOP 

27  OE  = QEO-r)QE 
SDO=0. 0 
CAL=0<’1  .E-3 

DO  3 IQE=I.NQE 
QE=QF  +OUE 
THETa=OE/57. 29578 


0''016900 

00017000 

00017100 

00017200 

00017300 

00017400 

00017500 

00017600 

00017700 

00017800 

00017900 

00018000 

00018100 

00018200 

00018300 

00018400 

00018500 

00018600 

00018700 

00018800 

00018900 

00019000 

00019100 

00019200 

00019300 

00019400 

00019500 

00019600 

00019700 

00019800 

00019900 

00020000 

00020100 

00020200 

00020300 

00020400 

00020500 

00020600 

00020700 

00020800 

00020900 

00021000 

00021100 

00021200 

00021300 

00021400 

00021500 

00021600 

00021700 

00021800 

00021900 

00022000 

00022100 
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TO=1.E10 

EM=ESO 

c ie^jo  is  a switch  signalling  end  of  Burning 

IEMP=0 

c assign  time  increment  for  integration. 
c 

dt=step 

C IWING  IS  A SWITCH  indicating  DEPLOYMENT  OF  WINGS 

1WING=0 

C NABLE  is  a switch  signaling  controls  DEPLOYED  FOR  GUIDED  FLIGHT, 
NABLt;  = 0 

c 

C GLIDE  - FUFO  - GLIDE  - FUFO  - GLIDE 

QiUKHi  calchlate  weight  constants  used  in  Smoothing  the  pitch  Rate 
CO***  USED  IN  G-BIAS  COMPUr at  ions 
SMLSTP=0. 1 

WTnUT 1=EXP (-SMLSTP/TMCCNG) 

wriNl=1.0-WTOUTl 

T0GH=TM0MMC+3.9 

AVThI)  = 0.0 

TM3=IMOMMC*3.0 

NTHSTP=61 

KTM  = fi 

C PK'INl  AtJU  LAHEL  run  de:script  lorj,  constants,  and 
c initial  conditions 

PR  IfiT  q,  r I TLE  , FFCTP  , VO  , EmO  , EMB  , D,  OE  » DT  , FC  , SPI  , V^^F  , ho  , HTERM,  VCW 

1,PWC,XWC 

9 FORMAT ( 1U1POAA/1h010XSHFFCTR13X2HV013Y2HM013X2HMB14X1HO/ 

1 IH  SFlS.b/lHO. 

2 6X9MUUAU  E'LtVrtXTHTM  STErRxGHTHRUSTSX  1 OHSP  I MPULSE9X6HV-W  I NO/ 

3 IH  SF  1 b . 6/ 1 HOAX  1 1 H IN  I T ALT  , F T4X  1 1 H H-fTM  ALT,FT15H  VEL  XW1ND,FT/S, 
A 12X,3hR^C,12X,3HXWC,/1H  bf-lS.C) 

PRINT  91 , VHXF ,VriYF 

91  F0PMaT(/8H  VitXF  = ,F12.2,20H  FT/SEC  VHYF  = ,F12,2,7H  FT/SEC, 

1/) 

PRINT  92, BELT, TM 

92  FORMaT(1H  19HTHPUST  RISE  TIME  = ,F12.4,5H  SEC,8H  TM  = ,F12.A, 

1 bH  SEC) 

C 

C 

C YO=lNniAL  ALTITUDE,  M 

C HO  = ALTirUI)E  READ  IN  IN  FT 

c ytekh=tepminal  altitude,  H 

c htfwm=terminal  altitude  read  in  in  Ft 

ir (lOPTY.NL. 1 ) go  to  29 
PR  ItiT  20  , XCG  , CLONG , AMOM  , BMOM 

26  format  ( IHO,  1 IHLOC  OF  CG  = , F 1 0 . A , 2X  , 3c,C  AL  , 1 AH  PROJ  LENGTH  =,El0.4, 

1 2x,3HCAL,lbH  AXIAL  M OF  I = , E 1 1 . b , £^X  , 7HKG  M**2, 

2 15H  TRANS  M OF  I = , E 1 1 . 5 , 2X  , 7HKG  M*<*2) 

29  IF(SPI.EQ.O.O)  GO  TO  60 

BRATr=FC/SPI 

IF (HRATE.EQ.0.0)  go  TO  60 
TB= (EMO-EMB) /URATE 


00022200 
00022300 
00022A00 
00022500 
00022600 
00022700 
00022800 
00022900 
00023000 
00023100 
00023200 
00023300 
00023AOO 
00023500 
00023600 
00023650 
00023700 
00023800 
00023900 
OOO2A0OO 
0002A020 
0002A0A0 
0002A060 
0002A1 00 
0002A200 
0002A300 
0002AA00 
0002Ab00 
0002A600 
0002A700 
0002A800 
0002A900 
00025000 
/0002bl00 
00025200 
00025300 
00025AOO 
00025500 
00025600 
00025700 
00025800 
00025900 
00026000 
00026100 
00026200 
00026300 
00026A00 
00026500 
00026600 
00026700 
00026800 
00026900 
00027000 
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initialize  time.  X.  X-DOT,  Y AND 

T = 0. 
x = o . 

Y = YO 
Z = 0. 

SPINrSPINO 
VHX= . 30^P*VHXF 
VHY=.30^I<<*VHYF 
XD=VFLO"COS(THFTa) ♦VHX 
YD  = VELO<-‘SIN(THLTa)  ♦VHY 
THFT  A = 57.29b78itATAN  (YD/XD) 

V = SQPT  (XLnt«2  + YI.Mn^2) 

ZD  = 0. 

SPIN(i  = 0.0 
XDO  = f!  .0 
YI)n  = 0.n 
zrin=n  .0 

YAW  = 0 . n 

alt=yu 

SPN(.=  p.O 

ISW=IBUPN(T,aLT,THETA.V) 

IF  ( I8W.EU. 1 ) GO  TO  70 
GO  TO  71 

70  TO=0.0 
dt=step/a. 

initialize  runge-kutta  subroutine 

71  call  RUNGLI  (U,V.P,p,2.FLIGHT) 

SOLVE  flight  eOUaTIONS  FOR  INITIAL 

call  FLIGhT(T,U.A) 

initialize  counter  for  determining 


0C027100 
00027200 
00027300 
00027A00 
00027500 
00027600 
00027700 
00027800 
00027900 
00028000 
00028100 
00028200 
00028300 
00028A00 
00028500 
00028600 
00028700 
00028800 
00028900 
00029000 
00029100 
00029200 
00029300 
00029A00 
00029500 
00029600 
00029700 
00029800 
00029900 
00030000 
00030100 
00030200 
00030300 
00030A00 
00030500 
00030600 
00030700 
00030800 
00030900 
00031000 
00031100 
00031200 
00031300 
00031A00 
00031500 
00031600 
00031700 
00031800 

Conditions,  ooo319oo 

00032000 

00032100 

00032200 

number  of  points  to  be  PLOTTED  00032300 


PRINT  AO.  TB 

AO  format ( 1H0.22HEFFECTIVE  BURN  TIME  = iFlO.A.AH  SEC) 
GO  TO  61 

60  TH  = 0'. 

BRaTE=0. 

COMPUTE  AUXILLIARY  CONSTANTS  AND  REDIMENSION  INPUTS 

61  CALS0=D*<>2<»1  .E-6 
0L0NG=D‘>1  .OE-3<»CLONG 
VELO  = 0.3UA8<»VO 

SUPVLL  IS  THE  SUPREMIUM  OF  PROJECTILE  VELOCITY. 
SUPVEL=VELO 

SUPAlT  is  SUPREMIUM  OF  PROJECTILE  ALTITUDE. 

supalt=o.o 
Vw=0.30A8«VwF 
VC»(  = 0.3046ovCW 
YO  = 0.30Ari<>HO 
YTFRM=0.3nAdOHTkRM 
RIERM=PE* YTEPM 

Y-UOT. 
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C AT  END  OF  TRAJECTORY  SOLUTION, 

C 

NRLOT=0 


initialize  counter  for  counting  lines  per  Page. 

LINF=0 


IPRiMTr-NPRINT 

YAWDEG=YAW<>DT0RAD 

GO  To  PRINT  OUT  INITIAL  CONDITIONS. 

GO  TO  4 

start  OF  solution  loop.  Save  last  values  of  flight  variables. 

5 XPrX 
RP  = R 
ZP  = Z 
XOP=XD 
YOP=YO 
VP  = V 

TfU  TcPrl  HETa 
CnACMP  = C.MACfl 
RFSISP=Pt5IS 
YAi^P  = YA'^'C'EG 
SPINP=SPIN 
SMAHGP=5MaRG 
STAF  AP  = STAFAC 
PSIP=PSI 
SR:.Gp  = SRnG 

Call  runge-kutta  subroutine  to  solve  flight  equations  from 

T TO  T+DT. 


42 

43 

44 


TO  Shorten  run  time  when  using  option  mmcsw,  change  time  step  to 

SMLSiP  FOR  NTMSTP  STEPS  ANU  THEN  RE^ur'E  INPUT  STEP  SIZE. 


IF 

NE. 

1)  GO 

TO 

‘*4 

IF  (T 

♦ STEP 

.LT 

• TM3) 

GO 

TO 

44 

IF  (T 

.LT.TM3) 

GO  TO 

4? 

KTM= 

kTM*  1 

IF  (K 

Tm.GE 

.NJT 

■vlSTP) 

GO 

TO 

43 

OT=SvlSTP 
GO  To  44 
0T  = T-'3-T 
GO  TO  44 
OT=STEP 

continue 

CALL  RUnGE2(T,DT) 


00''32400 

00032500 

00032600 

00032700 

00032800 

00032900 

00033000 

00033100 

00033200 

00033300 

00033400 

00033500 

00033600 

00033700 

00033600 

00033900 

00034000 

00034100 

00034200 

00034300 

00034400 

00034500 

00034600 

00034700 

00034800 

00034900 

00035000 

00035100 

00035200 

00035300 

00035400 

00035500 

00035600 

00035700 

00035800 

00035900 

00035902 

00035904 

00035906 

00035908 

00035910 

00035912 

00035914 

00035916 

00035918 

00035920 

00035922 

00035924 

00035926 

00036000 

00036100 

00036200 

00036300 
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C 


C 


c 

c 

c 


SAVE  POSITIONAL  COORDINATES  OF  PROJECTILE  FOR  LATER  PLOTTING  OF  00036400 


trajectory,  00036500 

00036600 

SAVE  Range  and  flight  time  in  arrays  for  subsequent  00036700 

POLYtjOMlAL  fit  00036800 

00036900 

SRNG=SQWT (X«X*(Y-Y0)**?*Z*Z)  00037000 

IF(T  .GE.TENABL)  IWING=I  00037100 

037200 

IF  (iMPLOT.FO.400)  go  TO  520  00037300 

IF  (SRNG.GT.sRNGEM)  go  to  520  00037400 

NPlOT=NPLOT ♦ I 00037500 

TS(NPL0T)=T  00037600 

RS  (Nt^LOT  ) =SRNG  00037700 

VS{Nt^L0T)=V  00037800 


*************00037900 


520  continue 

IF  ( 1GL11>E.E0. 1 ) SPIN  = CnEFL*DTORAD 
YAUOt  G = Y Aw*L)TOPAD 

IF ( IGLIUE.EO. 1 ) ST aF AC=YAW0EG+1HETA 
ThO= ( xn*YuD-YO*XDD) /V**2 
A V THI'  =W  T 1 N 1 * T HO  ♦ W T UU  T 1 * A V T riD 
IF  ( T .LT  .TOGB)  AVTriU  = 0.n 
IF  (V.GT.SUl^VEL)  SUPVFL  = V 
IF  (At  T.Gl  .SU^^AL^  ) SUPALT=ALT 
88  IF  ( ISa’.EO.O)  go  to  50 
IF  ( IF.NO.LiJ.  I ) GO  TO  51 
IF(T.&E.TO*TG*0ELT)  GO  TO  49 
GO  Tfi  51 

49  DT=STEP 
IEN'J=I 

CALL  BURN ( T , XMASS » thrust ) 

PRINT  80,  XMaSS.THRUST.ALT, V.T 
80  format ( IHO  ,9h  maSS  = ,F10.4,1IH  THRuST  = ,F10.4, 

I I3H  altitude  = ,FI0.2,I0H  SPEED  = ,F10.2,I0H  TIME  = ,Fl0.3) 
GO  TO  51 

50  lSw=IBURN ( T , alt, THETA  ,V) 

IF  (Ei'O.EU.EMH)  ISW  = 0 

IF  ( IS.'/.EO.O)  GO  TO  51 
TO=T 

PRINT  20,  TO 

20  FORMAT ( IHO , ISHBURN  STARTS  AT  ,FI0.4»4H  SEC) 

DT=STEP/4. 

51  IPRIrjT  = IPRlNT*l 

IF ( IPRINT.EO.O)  GO  TO  53 
GO  TO  52 

53  IPRINT=-NPR1NT 

ADVANCE  LINES  COUNTER  AND  CHECK  IF  Time  TO  EJECT  PAGE  AND  LABEL. 

4 L1NE=LINE*1 
DMY=TH0*UT0RAD 


00038000 

00038100 

00038200 

00038300 

00038400 

00038500 

00038600 

00038700 

00038800 

00038900 

00039000 

00039100 

00039200 

00039300 

00039400 

00039500 

00039600 

00039700 

00039800 

00039900 

00040000 

00040100 

00040200 

00040300 

00040400 

00040500 

00040600 

00040700 

00040800 

00040900 

00041000 

00041100 

00041200 

00041300 

00041400 

00041500 

00041600 
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Variable  is  the  angular 
output  of  choice 


116 

117 


c 

c 

c 

c 


FOP  COPPERHEAD  FLIGHTS  THE  DUMMY 

PITCH  rate  in  degrees/sec. 

OMY  IS  A DUMMY  variable  USE.D  FOR 
IF  (LiNt.LE.O)  GO  TO  6 
LINE=-50 

IF  ( IGLIDE.EO. 1 ) GO  TO  116 
WKITF  (6»7)  title 
GO  Tn  G 

WRITE  (6,117)  title  ^ , 

forma! ( IHl 20 AA/lH0,9HTIMEtSECS,7X3HX,M,5x.5HHAG,M«7X,3HZ,M,2X, 
1 8HX00T,M/S,2X8HYD0T , M/S , SXbHV , M/S  1 XdHM ACH  n0.2X7HDPAG,LB, 

•2  9H  theta, D,9H  PITCH, D,9H  DELTA, D, OH  800YA,D,9H  0THE,0/S) 

7 PORMaT(1h120AA/1H0,9HTIME,SECS,7X3HX,m5X5HHAG,M,7X,3HZ,M2X, 

1 HHXOOT ,M/52X6HY00T ,M/S5XSHV,H/S1X8BMaCH  NO . 2X7HDRAG , LB  * 

2 9H  theta, D,9H  YAW,U,9H  SPIN,R/^«9H  STA  FAC  ,8H  DUMMY  V) 

6 CONTINUE 

HAG=ALT-Y lERM 

PRINT  8,T , X,HaG,2,XD, YD,V,CMACH,RESIS,THETA, YAwDEG,SPIN,STAFAC 

8 forma! (IH  1F9.3,3F10. 1 ,3F 10. 1 ,F9.2,Eh 9.2,aF9.2) 

IF(T.GT.300.)  GO  TO  30 

RULE  FOR  stopping  SOLUTION  - STOP  WHLN  PROJECTILE  HITS  GROUND, 

52  IF ( .NOT . (R.Lf  .RTERM. ANO.THE! A.LT. 0. 0)  ) GO  TO  5 

INTERPOLAIE  SOLUTION  VAHlAdLES  FOR  R=RTERM 
YL  = Y TERM 

TE  = T-OT<»  (R-RTERM)  / (R-RP) 

Of L=(T-TE)/DT 
XE  = X-UFL‘'  (X-XP) 

ZE  = Z-DEL<»  (Z-ZP) 

XOE  = XU-DEL*  (XD-Xf)P) 

YDE=Y0-UEL* (YD-YDP) 

VE  = V-DEL<^  ( V-VP) 

THE  TAE  = THETA-nFL<^  (THETA-!  HE  TAP) 

CMACHE  = CMACh-DEL*>  (CMACH-CMACHP) 

RESISE=RESIS-0FL^ (RESIS-RESISP) 

YA‘WE  = YAV.'D'  G-OEL*^  ( YA wUEG- Y AWP ) 

SP I MK  =SP I N-UEL  * ( SP I N- SP I NP ) 
staf  ae=staf  AC-C»EL<MSTAF  AC-STaFAP) 

SmaRGE=5MARG-OEL* (SMAkG-SmaRGP) 

PS1E=PSI-DEL* (PSI-PSIP) 

SRNGE  =SRNG-DEL<>  (SRNG-SRNGP) 

YE=Yf -yTERM 
XPLOT (NPLOT) =XE 
YPLOT ( 1 ,NPL0T) =YE 

PRINT  OUT  SOLUTION  VARIABLES  FOR  Y=YTERM, 

PRINT  8,TE,XE,YE,ZE,XDE,YDE,VE,CMACHE,RESISE,THETAE,YAWE, 

1 SPINE, STaFaE, SRNGE 

SUPVFL=SUPVEL/0.30A8 
RANGE  = SURT  (XE«<»2  + ZE*<»2) 


000A1700 
000A1800 
000A1900 
000A2000 
000A21 00 
000A2200 
000A2300 
000^2400 
000^2500 
000^2600 
000A2700 
000A2800 
000^2900 
000A3000 
000A31 00 
000A3200 
000A3300 
,DMY000A3A00 
000A3500 
000A3600 
000A3700 
000A3800 
000A3900 
000A4000 
0004A  1 00 
00044200 
0004A300 
00044400 
00044500 
00044600 
00044700 
00044800 
00044900 
00045000 
00045100 
00045200 
00045300 
00045400 
00045500 
00045600 
00045700 
00045800 
00045900 
00046000 
00046100 
00046200 
00046300 
00046400 
00046500 
00046600 
00046700 
00046800 
00046900 
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RANGEsRANCZ" (1 .♦ ( R ANGE/RE ) **2/6 . ) / 1 a52 . 000A7000 

C PRINT  MAXIMUM  VELOCITY,  RANGE,  AND  ALTITUDE.  00047100 

PRINT  90,  SUPVEL,  RANGE, SUPALT  00047200 

90  format  ( 1H0,20HMAX  PROJ  velocity  = ,(115. 4, AH  F/S,  00047300 

1 3x,i2hmax  range  = ,fi5.a,iih  naut  Miles, 3x,iohmax  alt  = , ooo47aoo 

2 F15.A,8H  METERS)  00047500 

00047600 

Plot  the  trajectory  just  computed.-  00047700 

00047800 

Label  plot  with  title,  ue  and  vo.  00047900 

00048000 
00048100 

PRINT  1 0 , T I T L E , OE , VO  00048200 

10  format  ( 1 ho  1 5X?n A4/4H  OE  = F5.1,10H  DEG,  V0=F6.1,4H  F/S)  00048300 

3  CO^)TINUE  00048400 


call  polt it (R5,ts,nplot,3,o,cjrt,g,r, .true. ,SnEV, 

I 20HFLIGHT  TIME  VS  RANGE) 

call  polfit (rs.vs.nplot ,3,o,cjvt,h,himv, .true. ,SDEV, 

. 20HPROJ.  VEL.  VS.  RANGE) 

return  for  another  case. 

GO  Tu  101  00049200 

30  Call  exit  00049300 

^ 00049400 


00048600 

00048700 

00048800 

00048900 

00049000 
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SUnWOUTINE  SOUND ( A « G » RHO « V I SCO ) 

COMMON  EMOt£MB»SPl  « EC  » ATE  , DELT  t TO  « TR  , I SW  t V » THETA  * FFCTR  * C ALSO , 

1 VW,vCW,ALT,H, IEMO,CM-<CH,REYNLD,RESIS,cAL,nLONG,IOPTYtYAW,AMOM, 

2 HMOM»PSI  » WTAWF.  A « ISEP*NAULF  » IGLI  DE  » I w I NG * TEMABL  » A VTHD  * CDEEl  t 

3 MMCS'v  , tmommc 

common/snucom/caoens 

EUUIVAlENCE  {Y,ALT) 

subroutine  computes  the  speed  of  SOU\n  in  m/sec 

VERSUS  ALTITUDE  IN  METERS.  ALSO  COMPUTED  IS  THE 
ACCEl  EPATIOcj  due  TO  GPAVITy  IN  M/SEC/SEC  AND  THE 
AID  f,EFiSITY  IN  KG/M«<>3  AtoD  THE  ABSOLUTE  VISCOSITY 

OF  The  air  in  kg/h/sec.  note  that  Rlynold’S  number 
PER  meter  is  given  by  A*H-'HO<>EMACH/vibCO, 

*«**  REFERENCE:  BARNHART,  w' . TmE  STaNDARO  ATMOSPHERE  USED  BY  BRL 
<>«■•»*  FOR  trajectory  computations,  memo.  t<EPORT  NO.  1766.  JUNE  1966. 

OAIA  For  the  following  equations  wE'^E  extracted  from  the 
»»»*  U.S.  STANDARD  atmosphere.  1962. 

G = 9.P26<»  {6.37HE6/  (6.37RE6+Y)  ) **2 


0=6.356766E6+Y 

IF (Y.LE. 1 1019.07) 

GO 

TO 

1 

IF  (Y.LF.  .20063. 1?) 

GO 

TO 

2 

IF  (Y.LE. 32161.9) 

GO 

TO  3 

IF  (Y.LF.473b0.()9) 

GO 

TO 

4 

IF  (Y.LF_.  52428.88) 

GO 

TO 

b 

IF (Y.LE.' \591 .03) 

GO 

TO 

6 

IF  (Y.LE. 79994. 14) 

GO 

TO 

7 

RH0=0.a636«EXP (-0 

. 12?07E 

-3<>Y) 

TrTFi'iPFRAlURE  IN  DEGREES  KELVIN 
T=18n.65 

8 A=?0.053oS0RT (T) 

RHO  = RHO<>CAUEFiS 

VISCO  = O.OOA6  7«  (T  + IIO.  ) T/2 1 7 , 78  ) «<>  1 . b 

This  IS  the  SUTHERLAND  VISCOSITY  LA'^. 

return 

1 RhO=1 .22A999*Y<f  (-.  I 1 /60  3 3E-3*Y<»  ( . 4337  1 9E-8  + Y*  { - . 7A6  1 659E- 1 3 
1 ♦Y*  ( .S5376  0 3E-18-.95  72  7 27E-24<>Y)  ) ) ) 

T = ( 1 .831702ER-4. 1030H3E4<fY)/D 
GO  To  8 

2 RHO=l .R90 142* Y« (-.29401 14E-3*Y« ( . I 993974E-7 ♦ Y« ( - . 7637263E- 1 2 
1 ♦Y‘t(.  161S921E-I6-.  1476764E-?I<»Y)  ) ) ) 

T=216.65 
GO  To  8 

3 RHO=1.0156UY<M-.23b749E-3>Y*  ( . I 30807F-7  ♦ Y<>  { - . 38  1 965 1 E- 1 2 
1 ♦Y*(.579R729E-17-.3626654E-22<>Y)  ) ) ) 

T = ( I .^b00b8E9  + 6.b53416E3«Y) /D 
GO  10  8 

4 RHO=  1 . I 094A  + Y<*  {-.  1 I40029E-3  + Y'*>  ( ,481  74olE-8  + Y*>  {-,  I 039241E-12 
1 ♦Y<»(,lI38793E-I7-.b05213bE-23<*Y)  ) ) ) 

T=(8.839083E8+1.7938E40Y)/U 
GO  To  B 

5 RH0=,89  74979E-1*Y* (-.4179  05E-5  + Y* ( . 3529753E- I 0 ♦ Y« ( . 1 1 77 1 44E- 1 4 


00049500 

00049600 

00049700 

00049800 

00049900 

00050000 

00050100 

00050200 

00050300 

00050400 

00050b00 

00050600 

00050700 

00050800 

00050900 

00051000 

00051 100 

00051200 

00051300 

00051400 

00051500 

00051600 

00051700 

00051800 

00051900 

00052000 

00052100 

00052200 

00052300 

00052400 

00052500 

00052600 

00052700 

00052800 

00052900 

00053000 

00053100 

00053200 

00053300 

00053400 

00053500 

00053600 

00053700 

00053800 

00053900 

00054000 

00054100 

00054200 

00054300 

00054400 

00054b00 

00054600 

00054700 
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1 ♦Y*(-.2567072E-19+.UA9113E-24*Y)  ) ) ) 

T=270.65 
GO  TO  8 

6 RHO=,  10290R2t:-l  +7*  ( . 1 0«lB5JE-5  + Y*  (-•85236l9E-l0tY*  ( .2075003E-1A 
1 ♦Y*i-.2184B24E-19+.860425E-25*Y) ) ) ) 

T= (2.381 562E9-1 .?338h8E4*Y) /D 
GO  TO  8 

7 HHO=0.4636*EXP(-0.12?07E-3*Y) 

T=(3. 157088E9-2,h930^1EB*Y)/0 
GO  TO  8 ■ 

END 


0005A800 

00054900 

00055000 

00055100 

00055200 

00055300 

00055400 

00055500 

00055600 

00055700 

00055800 
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SUBROUTINE  FL I GHT ( T 1 ME , U . KUT T A ) 

DIMENSION  U(12) 

DIMENSION  T MACH (11)  «TKA(11)  »TK0YaW(11),TKL(11)  *TKM(11) « 

1 TkF ( 1 1 ) , TKT ( 1 1 ) , TKH ( 1 1 ) , TkS ( 1 1 ) » TCC  I 1) 

DIMENSION  TFMACH ( A ) , ALTRMM ( A ) * AAL  TRN ( A ) » TCN (4*  A ) • WTCN ( A) 
DIMENSION  TCAOl (A) , TCAD2 (A) , TDEL (4) »TaT(4,A)*WTAT(4) 

DAT  A GGA1N/0.50/,NTAAH/A. on/.UT5AH/b. 0 0/ 
data  DTORAD/S7. 29578/ 

DATA  RE/6.378E(S/,OMEGA/0. 729 1 5E-A/ , P 1 OFOR/ . 785398  1 /»  TWOG/ 1 9 . 584  1 
RErNDMlNAL  RADIUS  OF  THE  EARTH  AT  THE  EQUATOR  IN  METERS 
0MEGA=ANGULAR  velocity  OF  THE  EARTH  IN  RADIaNS/SEC 

table  of  equivalences 


1) 

= 

X 

2) 

= 

Y 

3) 

= 

Z 

4) 

SfMN 

b) 

= 

XOOT 

6) 

YDOT 

7) 

ZDOT 

8) 

SPINO 

9) 

= 

XIJBL 

10) 

1 — 

YUBL 

1 1 ) 

ZOHL 

12) 

1 — 

DUMMY 

c 

c 


EXTERNAL  CURAQ 

common  EhOiEMH  »SPI iFCfBRATE •DELT.TO* It3,ISW»V»THETA»FFCTR»CALSO» 
IVWO.VXWI alt iK, IEND.CmaCH,RE YMLD.RES1S,CAL.DL0NG, IOPTy* yaw, amom, 

2 BmOm,pSI *WTAREA,lSEPiNAtiLE»IGLlUE,IwlNG,TENABL*AVTHD*CDEFL» 

3 MmCSW,TM0MMC 

COMMON/COF  COM/XCG, SMaRG,EM, THID,PRNU  ,ALTRIM,CNATRM,GL1D£,EPSTHE 

1 , ST  AFAC  , yAWNU,  TKA  , TKUYAW  ,T><LiTKM,T^F,TKT,TKH,TKS,TCP,TMACH, 

2 NaRTBL*  TFMACH, ALTRMM,  AALTRM,TCN,  WTC|j,tCaD1 , TCAD2, TDEL*  tat ,WTAT 
COMMQN/W  1NC0M/RWC,XI'.C 

IF ( ISW.EO.O)  GO  TO  10 
IF  (TIMr.GT.TO*TP*DELT)  GO  TO  9 

Call  burn(Time,xmass, thrust) 

EM=XMASS 

Thrust  induced  drag 
ffc=ffctr<*thio 

H=A. AAa23»THRUST 
TERMX  = H<»U  (5) 

TERMY=H*U (6) 

12  VSO  = U(5)"<»2*U(G)*<*2*U(7)**2 
V=SURT (VSO) 

UP=RF*U(2) 

XSO  = U ( 1 ) <*<»2 
YSQ=UP«<*2 
ZSO  = U (3)  <»*2 
R=S0RT (XSO+YSQ+ZSO) 

ALT=R-RE 

ALT=U(2) ♦XS0/2,/RE 


00055900 

00058000 

00056100 

00056200 

00056300 

00056400 

00056500 

00056600 

8/00056700 

00056800 

00056900 

00057000 

00057100 

00057200 

00057300 

00057400 

00057500 

00057600 

00057700 

00057800 

00057900 

00058000 

00058100 

00058200 

00058300 

00058400 

00058500 

00058600 

00058700 

00058800 

00058900 

00059000 

00059100 

00059200 

00059300 

00059400 

00059500 

00059600 

00059700 

00059800 

00059900 

00060000 

00060100 

00060200 

00060300 

00060400 

00060500 

00060600 

00060700 

00060800 

00060900 

00061000 

00061100 


44 


oooo  ooo  o .o  oooo  on 


/ G LEVEL  21 


FLIGHT 


DATE  = 77143 


n/zz/z9 


DRCOSXsU ( 1 ) /R 

DRCOSY=UP/H 

DRC0S7=U(3)/R 

CALL • SOUND ( A*G»RH0» VISCO) 


THIS  GENERATES  SPEED  OF  SOUND » 6RA V I T Y , A I R DENSITY»  AND  VISCOSITY, 

IF(V.EO.n.O)  GO  TO  13 

TERMx=TERMX/V 

TERMY=TLRMY/V 

«<nt*COMPuTE  VELOCITY  OF  WIND  AS  A FUNCTION  OF  ALTITUDE, 
HARG=1.000<H)(2) 

VW=VWO+RwC*VWC (HaRG) 

VCW=VXW+XWC*VWC (HARG) 


vw=vwo 


vcw=vxw 

VRELSQ=  ((U(5)+VW)<»*2  + U(6)<*<*2*(U(7)  ♦VCW)  **2) 

VREL=S0RT ( VRELSQ) 


cmach=vre.l/a 

COr^PiiTE  REYNOLD'S  NUMBER  AND  SKIN  FRICTION  COEFFICIENT  (SFC), 
FRICT=0.0 

IT  (DLONG.LO.0.0)  GO  TO  4H 
REYNt,0  = ULON&'*A<»RHO"CMACH/VISCO 


ALR  = ALOG10 (Rf  YNLD) 


PwR  = O.Ob''ALR 

STc  = 0.4bS/ALR*'>2.bM/  ( 1 . ♦ 0 . 2*CMACH'*«^  ) **PWR 

IF  (ISEP.E-.'.O)  GO  TO  48 

FRICT=WTaREA<'SFC 

DEnOm=MASS  of  projectile  in  KG. 

48  DENOM  = n.4536'>EM 


DYNPRS  = 0.5'*RHO''VRELSCi 
IF ( lOPTY.EQ. 1 ) GO  TO  20 
YAW=0. 0 


YAWSO=0.0 

XKDYa«=0.0 

XTERmS=0.0 

YTERNS=0.0 

ZTERmS=0.0 


U(0)=0.0 

21  call  CDRAG(CMACH, drag. CAD) 

C0FnR6=FFC''DRAG*XKnYAw*  YAWSO  + FRICT  + CaD 

THIS  GENERATES  THE  CORRECTED  COEFFICIENT  OF  DRAG 


differential  EQUATIONS 
22  CONTINUE 

F0RM=PI0F  OR'^CALSU'^DYNPRS 
DR6=-C0FDRG*F0RM 


RESIS  IS  air  resistance  IN  POUNDS. 
RESIS=DRG  74.44823 

****  glide  - FUFO  - GLIDE  - FUFO  - GLIDE  - 
IF  ( IGLIDE.NE. 1 ) GO  TO  60 
C****  attitude  hold  logic  for  COPPERHEAD  AP77 


00061200 

00061300 

00061400 

00061500 

00061600 

00061700 

00061800 

00061900 

00062000 

00062100 

00062200 

00062300 

00062400 

00062500 

00062600 

00062700 

00062800 

00062900 

00063000 

00063100 

00063200 

00063300 

00063400 

00063500 

00063600 

00063700 

00063800 

00063900 

00064000 

00064100 

00064200 

00064300 

00064400 

00064500 

00064600 

00064700 

00064800 

00064900 

00065000 

00065100 

00065200 

00065300 

00065400 

00065500 

00065600 

00065700 

00065800 

00065900 

00066000 

00066100 

00066200 

00066300 

00066400 
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1200 

(;«■»*« 

C<nnn> 


1210 

Q-tt-O-ttO 

c 

C 


1212 

1220 


12U 


1230 

C 

C 


1200 

1200 


21 


THE=fiRSIN (U(6) /V) *DT0RaD 
IF (MMCSW.NE. 1 ) GO  TO  1260 
NOTE  THAT  OTHER  OPTIONS  FOR 


mechanizing  attitude-hold 


calculations  may  bE  EMPLOYED.  ThESL  aRE  FOUND  AFTER  LOC,  1260. 
toah=tmommc*dtoah 

TOAH  IS  The  time  at  which  the  attitude-hold  switch  is  thrown. 

IF (T IME.lt. TOAH)  go  to  60 
IF (NaSLE.EO. 1 ) GO  TO  1390 
T5Am=TM0MMC*UTS AH 

tsah  is  the  Time  at  which  the 
ARE  EXTENDED,  AND  THE  BODY  IS 
IF  (TlME.LT.TbAH)  GO  TO  1200 
NARLE=1 
IW1NG=1 
TEnaBL=T5aH 
GO  TO  1390 
CONTINUE 

calculation  of  CONTROL  DEFLECTION  fOr  G BIAS 
CDEFL=G(,Airj<»AHS  (AVIHD)  <>OTORAD 


GYRO  IS  freed,  wings 
commanded  to  FOLLOW  THE  GYRO, 


ARRAYS  TO  DETEI^MINE 
GO  TO  1230 


SET  UP  INTERPOLATION 
IF (CMACH.lt. TFMACH ( 1) ) 

DO  1210  1=2,4 

IF  (CmACH.LT.TFMACH'(  I)  ) GO  TO  1220 
CONT INUE 

LOCAL  MACH  NUMHFR  IS  BEYOND  THE  TABLt 
TR  IMMX  = AI.TRMM  (4) 

CADUTCADl  (4) 

CAn2=TCAn2 (4) 

DO  1212  J=1 ,4 


TRIM  attack  for  given 


00066500 

00066600 

00066700 

00066800 

00066900 

00067000 

00067100 

00067200 

00067300 

00067400 

00067500 

00067600 

00067700 

00067800 

00067900 

00068000 

00068100 

00068200 

00068300 

CON00068400 

00068500 

00068600 

00068700 

00068800 

00068900 

00069000 

00069100 

00069200 

00069300 


WTaT { J) =TaT (4, J) 

WTCN( J) =TCN(4,J) 

CONTINUE 
GO  To  1244 

FRaCT= (CMaCH-TFMaCH ( I-l ) ) / (TFMACH ( I ) -TFMaCH ( 1-1 ) ) 
DO  1214  J=l,4 

WTaT(J)=TaT(I-1,J)*FRACTo(7AT(I,J)-TaT(I-1,J)) 
wTCN(  J)  = TCN  ( I-l  , J)  ♦FRACH*  (TCN  ( I , J)  -T(.fj  ( I-l  , J)  ) 
COMT I NUE 

TRIt'MXrALTPMft  ( I-l  ) + FRACT<''  (ALTRMM  ( I ) "AL.TRMM  ( I-l  ) ) 
CaD1=TCA01  ( 1-1  ) ♦F pact TCADl ( I ) -TCADl  ( I- 1 ) ) 
CAD2=TCAD2 (I-l) ♦FRaCT<MTCAD2 ( I ) -TCaD2 ( I-l ) ) 

GO  TO  1244 
TRlMMX  = ALTRMf1  ( 1 ) 

CAL)1  = TCA01  ( 1 ) 

CA02=TCA02 ( 1 ) 

DO  1240  J=l,4 


WTAT ( J) =TaT ( 1 , J) 
WTCN(J)=TCN(1 ,J) 
CONTINUE 
CONTINUE 

BEGIM  INTERPOLATI 
DO  1246  J=2,4 
IF (CPEFL.lt. TDEL ( 


ON  FOR  TRIM  ATTACK 
J) ) GO  TO  1250 


00069400 
00069500 
00069600 
00069700 
00069800 
00069900 
00070000 
00070100 
00070200 
00070300 
00070400 
00070500 
00070600 
00070700 
00070800 
00070900 
00071000 
000711 00 
00071200 
00071300 
00071400 
00071500 
00071600 
00071700 
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1246 

1248 


CONTINUE 

CONTPOL  surface  deflection  required  EXCEEDS  CAPABILITY, 
WRITF  ERROR  message. 

write  (6tl248)  _ 

forma T ( IHO *' ERROR : CONTROL  DEFL . CALCULATED  FROM  G BIAS 


EXCEEDS 


ILIMIT* ) 

CALL  EXIT 
1250  COtriNUE 

FRACr2= (COEFL-TDEL  CJ-l ) ) / (TDEL ( J) -TDEL ( J-1  ) ) 

ALTRIM=WTAT  (J-1)  ♦Fr,ACT2<MWTAT  (J)-WTAT  (J-1)  ) 

calculate  the  normal  force  coeeficilnt 


DO  1252  J=2,4 

IF (ALTRIM.LT.AALTRM( J) ) GO  TO  1256 
1252  CONTINUE 

trim,  angle  Calculated  is  beyond  the  table. 

MESSAGE  AND  STOP. 

WRITE  (6»1254) 

1254  format ( IHO ♦' ERROR : TRIM  ANGLE  CALCULATED 


WRITE  ERROR 


S BEYOND  THE  TABLE' ) 


1256 


CALL  fXIT 
CONT  F.'UE 

FRACI 3= ( ALTRIM-AALTRM ( J-1 ) ) / ( AALTRM ( J) -AALTRM ( J-1 ) ) 

C^l=V^TCN  ( J-1  ) •fFRACT3«  (WTCN  ( J)  - WTCM  ( J-  1 ) ) 

UPDATE  INITIAL  BODY  ATTITUDE 

ADf'OYO  = THE  + aLTR  IM 

GO  Calculate  aero  forces 

Go  TO  1500 
CO, NT  INUt 

IF  (TIME.LT.TENA&L) 

IF (NaHlE.NE. 1 ) GO 
GO  TO  1390 
CONTINUE 
TIME  .GE.  TENARL 
IF (NaBLE.EO. 1 ) GO 
GPEPS^GLlDE  + e PSTHE 
IF (The .GT .GPEPS)  GO  TO  60 

CALCoLaTION  for  INITIAL  CONDITIONS  EON  ATTITUDE  HOLD 

attitude  hold  starts  with  enable  SEI  T()  1,  ENTER  ONLY  ONCE  FOR 


1260 


66 


GO  TO  66 
TO  1300 


TO  1390 


1300  NAPLf  =l 

FIRSl  iterate  FOR  CNO  AND  ALTRIM 
KOUNTl=l 

CNO=DENOH<»TWOG»U  (5)  /V/2.  O/EORM 

This  is  The  required  normal  force  coef . to  preserve  glide  angle. 

COmPUTL  max  trim  At'lGLE  , 

IF (CmaCH.lt. TFMACH( 1 ) ) GO  TO  1330 
DO  1310  1=2*4 

IF(CmACH.LT.TFMACH(1) ) 60  TO  1320 

1310  continue 

(;««««  local  MACH  NUMBER  IS  BEYOND  THE  TABLE 
TR  lMtiX  = ALTRMM  (4) 

CAD1=TCAU1 (4) 

CA02  = TCAL)2  (4) 

DO  1312  J=1 *4 
WTCN( J) =TCN(4,J) 


00071800 

00071900 

00072000 

00072100 

00072200 

00072300 

00072400 

00072500 

00072600 

00072700 

00072800 

00072900 

00073000 

00073100 

00073200 

00073300 

00073400 

00073500 

00073600 

00073700 

00073800 

00073900 

00074000 

00074100 

00074200 

00074300 

00074400 

00074500 

00074600 

00074700 

00074800 

00074900 

00075000 

00075100 

00075200 

00075300 

IC00075400 

00075500 

00075600 

00075700 

00075800 

00075900 

00076000 

00076100 

00076200 

00076300 

00076400 

00076500 

00076600 

00076700 

00076800 

00076900 

00077000 
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1312 

1320 


1314 


1330 

C 

C 


1340 

1350 


1360 


1366 


1370 


1374 


1376 

1377 

1378 


WTaT ( J) =VAT (4, J) 

continue 

GO  TO  1350 

FHaCT=(CMACH-TFNACH(I-1) )/(TFMACH(I)-TFMACH(I-1)  ) 
DO  1314  J=l*4 

WTCN( J) =rCN( 1-1 , J) ♦F4ACT* (TCN (I,J)-TCN(I-1»J) ) 
WTAT(J)=TaT(1-1iJ)*FWACT*(TAT(1»J)-IaT (I-1»J) ) 

continue 

TR1N‘-;X  = alTHmm  (1-1)*FHACT«(aLTRHM(1)~ALTKmm(I-1)) 
CADI  =TC ADI  (I-l)  ♦F«ACT<»  (TCADl  (1)-TCAD1  ( 1-1)  ) 

CA02  = TCAD2  ( 1-1)  ♦FkACT<»  ( TCA02  ( I ) -TCAU2  ( I-l  ) ) 

GO  TO  1350 
THlt^MX  = ALTRMM  ( 1 ) 

CAOl  = TCADl ( 1 ) 

CA02=TCAr)2  ( 1 ) 

DO  1340  J=1.4 
WTCN ( J) =TCN ( 1 ♦ J) 

'a'TAT  ( J)  = TAT  ( 1 , J) 

COnT INUt 

continue 

IMEKPOLATE  IN  v.TCN(J)  FOR  THE 
ANO  CONTROL  SURFACE  DEFLECTION 
DO  1360  J=2^4 

GO  TO  1370 


instantaneous 
at  trim. 


TRIM*  ALTRIM 


IF  (C'lO.LT  .WTCN  ( J)  ) 

CONTINUE 

REOUIRED  NORMAL  FORCE  EXCEEDS  CAPABILITY.  PRINT  NOTICE  AND  SET 
ALTRIM  TO  TO  MAXIMUM. 

ALTKIM=TR IMMX 
WRITt  (6.1366) 

FORI'aT  ( IHO.  ' Ir^ilT  lAL  normal  FORCE  REQUIRED  FOR  DESIRED  GLIDE 
s capability • ) 

REPLACE  CNO  .vITH  MAX  REALIZABLE 

CNO=WTCN (3) ♦ ( TCN (4) -WTCN ( 3) ) * ( T R 1 mMX- A ALTRM ( 3 ) ) / 

(AAl  TRM(4)-AALTRM(3)) 

GO  TO  1374 
COr-'T  INUE 

INTERPOLATE  FOR  INITIAL  TRIM  ANGLE  REOUIREO  FOR  NORMAL  FORCE. 

F R ACT  2= (CNO- WTCN (J-1)  )/(WTCN(J)-wTCN(J-l)  ) 

ALTRIM=AALTRM  ( J-1  ) ♦FRACT2*^  ( AALTRM  ( J)  - AALTRM  ( J-1  ) ) 

CONT INUE 

INTERPOLATE  FOR  CONTROL  SURFACE  DEFLECTION 
on  1376  J=2.4 

1F(ALTR1M.LT.WTaT ( J) ) GO  TO  1378 
COI.TINUE 

ALTR'IM  is  oEYOND  the  table.  write  E,-P0R  message  AND  STOP, 

WRITt  (6.1377) 

FORMAT ( IHO. 'ERROR.  ALTRIM  BEYOND  THE  TABLE.  CHECK  INPUT 
CALL  EXIT 

continue 

FRACT3=( ALTRIM-WTAT (J-1) )/(wTAT(J)-WTaT (J-1) ) 

CDEF  L = TDEL  ( J-  1 ) ♦FRACT3''  ( TOLL  ( J ) - T DEL  ( J- 1 ) ) 

K0UNT1=K0UNT1*1 
IF (KoUnTI .GT.2)  GO  TO  1375 


00077100 
00077200 
00077300 
00077400 
00077500 
00077600 
00077700 
00077800 
00077900 
00078000 
00078100 
00078200 
00078300 
00078400 
00078500 
00078600 
00078700 
00078800 
00078900 
00079000 
00079100 
00079200 
00079300 
00079400 
00079500 
00079600 
00079700 
00079800 
00079900 
EXCEEDOOOHOOOO 
00080100 
00080200 
00080300 
00080400 
00080500 
00080600 
00080700 
OOiiBOBOO 
00080900 
00081000 
00081100 
00081200 
00081300 
00081400 
00081500 
00081600 
PARAMS’ ) 00081700 
00081800 


00081900 

00082000 

00082100 

00082200 

00082300 
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CNO=CNO/COS (ALTRIM/DTORAD) 

GO  TO  1350 

EfiO  OF  CALCULATION!  OF  INITIAL  TRIM  ANGLE  FOR  GLIDE  OPTION  1, 

NO.N  calculate  initial  body  attitude* 

1380  AH0()Y0  = THE  + ALTRI'-! 


CN=C\Q 
GO  To  1500 

C««<n>  SU‘j5F,0UEM1  calculation  of  ALTRIM  FOf<  attitude  hold 

At  TRlM  = AciODYO-THE 

IF  (Ci'ACH.LT.TFMACHd  ) ) GO  TO  1430 
DO  1410  1=2*4 

IF(CfiACh.LT.TFMACH(I)  ) GO  TO  1420 
1410  continue 

(;««««  local  MACh  NUMBER  BEYOND  THE  TAULE 

TRit'tiX  = ALTRMM(4) 

C CA01=TCAD1 (4) 

C CAD2=TCAi)2  (4) 

DO  1412  J= 1*4 
wTCN  ( J)  = TC!J  (4*  J) 

WTaT(J)=TAT (4, J) 

14  12  COiiTlNUL 


GO  10  1450 
1420  F RaCT= (CmaCH- 
DO  I4l4  J=1 *4 
WTCN  (v!)  =TCN  ( I 
WTaT ( J) =TAT ( I 
1414  CONTINUE 

TRIMmX=ALTRMM 
C CAD1=TCAU1 ( I- 

C CAD2=TCAD2 ( I- 

GO  TO  1450 
1430  TR  Im*-X  = ALTRMM 
CAD1=TCA01 ( 1 ) 
CAD2=TCAU2 ( 1 ) 
DO  1440  J=l*4 
WTCN ( J) =TCN ( 1 

wtat ( J) =1  at ( 1 
CONT INUE 
CONI INUE 
INTE'^POLATE  I 
DO  1460  J=2*4 
IF  (ALlRlM.LT. 
COfjTlNUE 
REOUlRFU  TRIM 
BE  LIMITED  TO 
altrim=trimmx 
CN=WTCN(3)  + ('/ 
1 (aaltrm(-)-a 
GO  TO  1474 
1470  continue 

c««««  interpolate  F 


1440 

1450 


1^*60 


TFMaCH ( I - 1 ) ) / ( TFMACH ( I ) -TFMACH ( I - 1 ) ) 

-1  * J) +FRACT4 (TCN ( I * J) -TCN ( I-l * J) ) 

-1  * J) *FRACT4 (TAT (I,J)-!aT(I-1*J)) 

(I-l) ♦FRACT*(ALTRMM(I)-ALTRMM(I-1) ) 

1 ) + FRACr<‘  (TCAOI  ( I ) -TCaUI  ( I-l  ) ) 
1)+FRACT4(TCAD2(I)-TCAU2(I-1)  ) 


(1) 


* J) 
f J ) 


N *TCN(J)  FOR  THE  NORMAL  FORCE  COEF, 
AALTRM( J) ) GO  TO  1470 

angle  exceeds  max  tabulated  value. 

TRIMMX. 

TCN (4) -WTCN (3) ) « ( ALTRIM-AALTRM (3) ) / 
ALTRM (3) ) 


OR  CN  within  the  Table 


CN. 


TRIM  WILL 


00082400 

00082500 

00082600 

00082700 

00082800 

00082900 

00083000 

00083100 

00083200 

00083300 

00083400 

00083500 

00083600 

00083700 

00083800 

00083900 

00084000 

00084100 

00084200 

00084300 

00084400 

00084500 

00084600 

00084700 

00084800 

00084900 

00085000 

00085100 

00085200 

00085300 

00085400 

00065500 

00085600 

00085700 

00085800 

00085900 

00086000 

00'86100 

00086200 

00086300 

00086400 

00086500 

00086600 

00086700 

00086800 

00086900 

00087000 

00087100 

00087200 

00087300 

00087400 

00087500 

00087600 


. ^9 

9 


on  non  on 


G LEVEL  21 


FLIGHT 


DATE  = 771A3 


17/22/29 


FRACT2= (ALTRIM-AALTRM ( J-1) ) / ( AALTHM ( J) -AALTRM ( J-1 ) ) 

CN=WTCM(J-1) ♦FRACT2*(WTCN(J)-wTCN(J-1) ) 

1474  COfjTiHUE 

interpolate  for  control  surface  deflection 

DO  1476  J=2»4 

IF(ALTRlM.LT.'rtTAT ( J) ) GO  TO  147B 
1476  continue 

c«<>*«  altkim  is  beyond  the  Table,  write  error  message  and  stop. 

WRITE  (6.1377) 

Call  exit 

1478  CONTINUE 

FRACT3=  ( ALTRIM-V;  IAT(J-1))/(WTAT(J)-WTaT(J-1)) 

COEFl.  = TUEL  ( J-1 ) ♦FRACI3<»(TOFL(J)-TDEL(j-l)  ) 

end  of  normal  force  coefficient  calculations 
program  reentry  for  ATTITUUE-HOLO  CALCULATIONS 
isoo  continue 

CONVERT  angle  TO  RADIANS 
C0EFL  = CDEFL/L'T0RA0 
ALTRIM=ALTRlH/nTORAO 

C<nntu  PLACE  TRIM  aNGLI  IN  YAW  POSITION  FOR  PRINTOUT. 

YAW=ALTRIM 
FN=CN<^FURM 
SRJAI.=S1N  ( ALTRIt-i) 

COSAL=COS ( ALTRIM) 

Calculate  axial  force  increment  produced  by  fin  deflection  at  t 

UELCa=CAD1  ( 0 .4<»AL  TRIM-CDEFL  ) **2-CAD£.’o’.lTRIM<^*3 
DFLCA=COEFI/6.0 
FDI2=nELCA*FURM 
FL=FN<»C0Sal-FDI2"SINAL 
F0I=-EN<^SINAL-F0I2<^C0SAL 
sinth=u(G)/v 
cosrH=u (5) /V 

TEPMX  = TLRMX-FL<*SlNTh  + FDI<^COST  H 
TERMY  = TERMY  + FL*C0STH  + FDI<>S1NTH 

60  COr.TiNUE 

IF (NauLE.EQ. 1 ) IWING=1 
END  OF  ATTITUDE-HOLO  CALCULATIONS 

U(  10)  = (DPG«U  (6)  / VREL  + TERMY)  /DENOM-G*URCOSY  + 0.b3166E-8<>UP 
1 ■♦?.*C)MLGa<*U  (7)  ♦YTERMS 

U(9)  = (ORG<^  (U(5)  ♦V.v)  /VREL  ♦TEPMX)  /DENO^'-G«DRCOSX♦XTERMS 
U(  1 1 ) =-2."0MEGA"U (6) -G*DRC0S7«ZTERm8  + DRG<>  (U(7) ♦ vCW) /VREL/DENOM 
U( 12) =0.0 

AC  IS  THE  ACCELERATION  OF  THE  PROJECTILE  ALONG  THE  TRAJECTORY. 

AC=  (U(5)*U(y)+U(6)<*U(10)+U(7)<>U(ll))/V 
lEdOPTY.EO.O)  GO  TO  14 
U (B)  =-RH0<*CALS0<n»2/AM0n<*XKA<»U  ( 4)  <»V 
14  IE  (KUTTA.EO.4)  THETA=ARSIN(U(6)/V)"57.29578 

return 

13  U(9)=0. 

U(10)=-G 

THLTA=90. 


00087700 

00087800 

00087900 

00088000 

00088100 

00088200 

00088300 

00088400 

00088500 

00088600 

00088700 

00088800 

00088900 

00089000 

00089100 

00089200 

00089300 

00089400 

00089500 

00089600 

00089700 

00089800 

00089900 

00090000 

RIM00090100 

00090200 

00090300 

00090400 

00090500 

00090600 

00090700 

00090800 

00090900 

00091000 

00091100 

00091200 

00091300 

00091400 

00091500 

00091600 

00091700 

00091800 

00091900 

00092000 

00092100 

00092200 

00092300 

00092400 

00092500 

00092600 

00092700 

00092800 

00092900 
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RETURN 
9 EM=EMH 
GO  TO  11 

10  EMsEMO 

11  TERMX=o. 

TERMY=o . 
ffc=ffctr 
on  To  12 

20  CALL  ACOEFS (CMACH » Y A A « XK A » XKDYAW » XKL » XKM » XKF * XKT » 
IXKHfXKStXCP) 

VXRL=U (5) ♦VW 
VZRL=U (7) ♦VCW 

VRLS0  = VXRL<><>2*IJ  (6)  «o2*VZRL«*2 
VRL=S0RT ( VRLSQ) 

C TEST  for  dynamic  SIADILITY. 

HOTTOM=8.n<>!5MOM«DYI]Pr<'S«CAL<»*3oXKM 
STaF aC=  (U(4)  ■»AM0M)  <>«2/HOTTOM 
If  (STAFAC.LE. 1 .0)  GO  TO  25 
COMPUTE'  the  yawing  FE'EOUENCY 

yAWNn  = AMOM/hMOMOS.jRT  ( 1 . - 1 . /ST  AF  AC  ) <^U  ( A ) /6 . 2832 
C«<n>«  COMI’UTE  THE  iiVEPTURNING  mOI’ENT  AND  RRECESSIONAL  FREQUENCY 
OT  fv'-'Ovi  = 2 . <>Xkm0CAL“<>3<M)Y  OPRS 
PWiJU  = OTNiiOM/AMOfi/U  (A)  /6.2832 
C CO'-PUTf  YAW  OF  REPOSt. 

, ALPHA-i  = RMOi>'CAL'‘WPLSiJ<>  ( X KL'’’ XK M* VPL SO ♦ C ALSQ'*>XKF*XK T*U  (4)  **2) 

IF  ( ArS (ALPHAM)  .LT. 1 .h-20)  GO  TO  25 
ALPhaa  = AMO:-i<-XKL.'“U  (4)  /CalSn/ALPHAB 
ALPHAH=DENOmvXKT«U (A) /ALPHAB 
AME.  = ALPHAb-ALPHAA 

ALPHAX  = AMr3o  (U  (G)  ( 1 1 ) -VZR'LOU  (10))  -ALPHAP* VZRL*G 

ALPU''' Y = AMd<>  (VZRL0U(9)-VXK'L<^U(11)  ) 

ALf'r3AZ  = AMH<''(VXPL^HJ(10)-U(b)*U(9))+ALPhAH''>VXRL*G 
YAWSN=ALPhAX*«2+ alpha Y<h>2+ALPHAZoo2 
YAW  = SQRT  ( YAftSCJ) 

IF (YAW. GT. 1 .5708)  GO  TO  25 

ARG=  ( VXRL<»AI.PHAZ-V/kL<>ALPHAX  ) OVRL 

A(7(S1  = (VXRL*ALPHAY-U(G)oalP|(AX)«’VXPL“(U(6)<*ALPHAZ-VZRL*ALPHAY) 
IF  ( APS ( aRGI ) .LF  . 1 . 0E-2n ) GO  TO  50 
PS1=57.3oaTaN(ARG/ARG1 ) 

GO  Tc  53 

50  IF  ( Af.  (i'»AKGl  ) 51,52,52 

51  PSl=-yo. 

GO  TO  53 

52  PSI=P0. 

53  continue 

PSI=0R1ENTATI0N  OF  YAW.  This  is  the  angle  between  the  plane 
containing  BOTH  THE  VELOCITY  AND  YAW  VECTORS  AND  A VERTICAL 
plane  CONTAiriir-'G  THE  VELOCITY  VECTOR.  IT  IS  MEASURED 
CLOCKWISE  FROM  THE  VERTICAL  PLANE. 


END  OF  COMPUTATION 
DKFN=CaL*XKFou (4) 


OF  YAW. 


00093000 

00093100 

00093200 

00093300 

00093400 

00093500 

00093600 

00093700 

00093800 

00093900 

00094000 

00094100 

00094200 

00094300 

00094400 

00094500 

00094600 

00094700 

00094800 

00094900 

00095000 

00095100 

00095200 

00095300 

00095400 

00095500 

00O95600 

00095700 

00095800 

00095900 

00096000 

00096100 

00096200 

00096300 

00096400 

00096500 

«VZRL00096600 

000'''')700 

00096800 

00096900 

00097000 

00097100 

00097200 

00097300 

00097400 

00097500 

00097600 

00097700 

00097800 

00097900 

00098000 

00098100 
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XKLVSQ=XKL*VRLSO 
HUSO  = RHO*CaLS(Vi)ENOM 

XTFWmS=ROSOo  (XKLV50<*ALPHAX*DKFN*  ( ALPriAY«VZPL-ALPHAZ*U  (6)  ) ) 
YTFP''’S=HObOo  ( XKLVSO*  alpha  Y*DKFN»  ( ALPhAZ*VXRL-ALPHAX*VZRL)  ) 
ZTEHmS  = R()SQ»  (XKLVSQ*ALPHAZ*DKFN»  (ALPHaX‘»U  (6)  -ALPHAY*VXRL)  ) 

XKOYAW  = 2,5At)'47*XKDYA'w 

IFCTAVJ.LT.O.eP)  GO  TO  21 
PRINT  55»RH0»XTtRKS» YTFRMS,ZTERMS 
55  FORMAT (IH  .1P4F10.S) 

GO  TO  21 

25  PRINT  ?6,STAFAC 

26  format  (1H0»40HUNSTABLE  PROJECTILE  S>TaBILITY  FACTOR  = iFlO.4) 

call  exit 

END 


00C98300 

00098400 

00098500 

00098600 

00098700 

00098800 

00098900 

00099000 

00099100 

00099200 

00099300 

00099400 

00099500 

00099600 

00099700 

00099800 
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SUBROUTINE  CDRaG (EMaCH.DRAG»CAD) 

PROGRAM  COMPUTES  THE  COEFFICIENT  OF  DRAG  VERSUS  MACH  NUMBER 
AND  THE  COEFFICIENT  INCREMENT  DUE  TU  WINGS. 

DIMENSION  XMTBL(12) .cnT8L(12) .C0T8L<12) 

DATA  COOKD/2. 546479/  ^ 

COMMON  EMOiEmb.SPI.i-C.BRATE.DELT.TO’TB.ISW.V.THETA.FFCTRiCALSO. 

1 VW.VCWiALT.t-’,  IEND.CMACH,RF,YNLD«RESiS»CAL»DLONG.  10PTY.YAW.AM0M. 

2 BMOm.PSI  . wT  ARF  a . I SF  P.MAbLE.  IGLIOE  . I ING.  TENABL  . AVTHD.  CUEFL  . 

3 MMCSW.TMOMMC 

C0''1M0N/DkGC0M/  xmtbl.cdtbl.cotbl  .ntbl 

CAI)  = 0. 0 

■ If  (MTBL.NE.O)  GO  TO  5 
IF  (E'‘'ACH.LF..O.PO)  GO  TO  1 
IF  (E^'ACH.LE.  1 . 10)  GO  TO  3 
If  (LMACH.Lt.3.0)  GO  TO  4 
EM3=FMaCH-3. 0 

DWaG=0 . 09  + Lm3<M  - 0 . 02*  0 . 002<‘EM3) 

0RaG  = CDOKD<>DRAG 
RLlUPN 

1 ()PaG  = 0.0589 

DRaG  = CDOKU<>URAG 
R t 1 U P ^4 

3 C=10.« (LHACH-n.8) 

DRAG  = 0.07  7 36'-C<f<>3<»EXP  (-C)  +0.0589 
URaG=COORD<>DRAG 

return 

4 DRaG  = 0.2I547  + EMACH<^  (-0.05134+0. 00317«EMACH) 

ORAG  = CDOKD<>ORAG 

RETURN 

5 DO  6 J=l.NTfjL 

IF  (t haCH.lt. XMTBL (J) ) GO  TO  8 

6 CONTINUE 
ti  JL=J-1 

IF ( Jl  .LC.O)  GO  TO  12 

f RAC=  (FMaCH-XMTGL  ( JL)  ) / (XMTBL  (J)  -XM.TdL  'JD  ) 

CO=COThL(JL)  + (COTr(L(J)-CDThL(  JL)  )*f  « AC 
CAO=0.0 

IF ( IWInG.NE. 1 ) GO  TO  lO 

CaO=COTUL  ( JL  ) + (COTliL  ( J)  -COTBL  ( JL)  ) *FRaC 

10  continue 
r)RAG=cn 
cad=cao 

RETURN 
12  CONTINUE 

ORAG  = CnTriL  ( 1 ) 

CAD  = 0. 

IF  (IWlNG.EO.l)  CAD  = COTBL(l) 

RETURN 

END 


0C099900 
00100000 
00100100 
00100200 
00100300 
00100400 
00100500 
00100600 
00100700 
001 00800 
00100900 
00101000 
00101100 
00101200 
00101300 
00101400 
00101500 
00101600 
00101700 
00101800 
00101900 
00102000 
00102100 
00102200 
00102300 
00102400 
00102500 
00102600 
00102700 
00102800 
00102900 
00103000 
00103100 
00103200 
00103300 
00103400 
00103500 
00103600 
00103700 
00103800 
00103900 
00104000 
00104100 
00104200 
00104300 
00104400 
00104420 
00104450 
00104500 
00104700 
00104800 


53 


i 


G LEVEL  21 


ACOEFS 


DATE  = 77143 


11/22/29 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUpROUT INF  ACOEF5(EMACH»YAw»XKAfXKDYAWfXKLfXKM»XKF»XKTfXKH»XKS»XCP00l'‘^4900 

‘ 00105000 

D I^iEMSlON  THACH(ll)  fTKA(ll)  »TKOYAW(ll)  fTKL(ll)  ♦TKM(ll)  ♦TKF(ll)  f 
1 iKTdl)  fTKH(ll)  ,TK5(11)  •TCP(ll) 

DIMENSION  TCAUl  (4)  t TCAD2 (4)  » TUEL (4)  f T AT  14,4)  , WTaT (4) 

COMMON/ COFCOf-VXCG,SMAK('i»  EM,  THIOt  PWmD  , ALT  R I M , CN  A TRM  f GL  1 DE  , EPSTHE 

1 ,ST/'FACfYAWriU,TKA*TKDYA'*fTKLtTKri»T^F  fTKT  fTKHtTKS»TCP»TMACH» 

2 NAK  FbLf  TFMACH,  ALTRM^if  AALTRM,TCNf  WTCN,TCAD1  »TCaD2»TDEL»TAT»WTAT 
EMaCH  = MACH  NUMhFR 

XKA  = SPIN  DAMPING  MOMENT  COEFFICIENT 
XKDYAW  = YAW  DRAG  COEFFICIENT 


XKL 

XKM 

XKF 

XKT 

XKH 

XKS 

XCP 


lift  force  coeeeicilnt 

OVER turn TNG  MOMENT  COEEF IClENl 
MAGNUS  FORCE  COEFFICIENT 
MAGNUS  MOMENT  COEFFICIENT 
DAMPING  moment  COEFrlCIENT 

pitching  force  coefficient 

CENTER  OF  PRESSURE  AFT  OF  NOSF  IN  CALIPERS 


FOR  OEPENOFMCE  OF  ACOEFS  UPON  YAW  St.E  PPL  MEMO.  RPT*  NO*  2023 
RELAIIVL  to  T3R7  TYPE  PROJtClILE. 

XKT=-0. 14^0.0576<'^  ( EMACH- 1 . 25  ) <^<^2 

XKT=o.n 
XKF  =(>.  157 

SYAV.=SP4  { YAW)  ^^2 

GO  TO  50 

51  CONTINUE 

00  60  J=1^^JA^^TBL 
IF  (t*'AC»ULT  .TMACH  ( J)  ) GO  TO  70 
60  COMlINUE 
70  JL=J-1 

FI<'AC=  (EMACH-TMACH  ( JL)  ) / ( TMaCH  ( J)  -TmACH  ( JL)  ) 

IF  (Tk  A ( J)  .EQ.0.0)  GO  TO  b2 

XKA=TKA  ( JL)  ♦ (TKA  ( J)  -TKA  ( JL)  ) *FT?AC 

52  IF  (ThOYAiv  ( J)  .LO.0.0)  GO  TO  53 

XKOYa»“  = TKDYAw(JL)*(TKL)YAw(J)-Ti<UYAW(JL))*FRAC 

53  IF  ( If L (J)  .EC.  0.0)  GO  TO  54 

XKL  = TKL  ( JL)  ♦ ( TKI.  ( J)  - IKL  ( JL)  ) *FHAC 

54  IF  (TKfM  J)  .FO.O.O)  GO  TO  55 

XF,M  = TKM  ( JL)  ♦ (TKM  ( J)  -TKM  ( JL)  ) «FPAC 

55  IF  (TKF ( J)  .EO.O.O)  GO  TO  56 
XKF=TKF (JL)*(TKF (J)-TKF(JL))*FRAC 

56  IF ( IKT ( J) .LO.O.O)  GO  TO  57 
XKI=TKT  (JL)*  (TKT  (J)-TKT  (JL)  )<fFRAC 

57  IF  (TfH(J)  .EO.O.O)  GO  TO  58 
XKH=IKH ( JL) ♦ ( TKH ( J ) - TKH ( JL ) ) *F R AC 

58  IF  (TkS( J)  .EQ.O.O)  GO  TO  59 
XKS=TKS ( JL) ♦ ( TKS ( J) -TkS ( JL) ) *FRAC 

59  IF  (Tkm(J)  .EQ.O.O)  GO  TO  62 
IF (XKL. EO.O.O)  CALL  EXIT 
SMaRG=-XKM/XKL 
XCP=XCG*SMARG 

RETURN 


00105100 
00105200 
00105300 
00105400 
00105500 
00105600 
00105700 
00105800 
00105900 
001 06000 
00106100 
00106200 
00106300 
00106400 
00106500 
00106600 
001 06700 
00106800 
00106900 
001 07000 
00107100 
00107200 
00107300 
00107400 
00107500 
00107600 
00107700 
00107800 
00107900 
OOlOROOO 
00108100 
00108200 
00108300 
00108400 
00108500 
00100600 
00108700 

ooloeeoo 

00108900 
00109000 
00109100 
00109200 
00109300 
00109400 
00109500 
00109600 
00109700 
00109800 
00109900 
001  10000 
001 1 01 00 
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62  IF(TCP(J) .EQ.0.0)  GO  TO  63 

XCP  = TCP  ( JL)  ♦ (TCP  ( J)  -TCP  ( JL)  ) <^FRAC 

63  SMflPCi  = XCP-XCG 
XKM  = -XKL<>SMAPG 
PETUPN 

50  CONTINUE 

IF (FNACH.LE.O.R)  go  TO  10 
IF  (E''.ACH.LE.0.9)  GO  TO  20 
IF (LOACH. LE. 1 .0)  GO  TO  30 
IF  (E*' ACH.LE.  1 . 1 ) GO  TO  35 
IF (EhACH.LF. 1 .30)  GO  To  AO 
IF  (f>’ACH.GT.  1 .5)  GO  TO  A5 
C VALID  FUR  EHACH  GTR  O.H  AND  LT  1,5 

5 XKA  = 0.0  03B*().  0 02"EXP  (-1 .5<>  (EMACH-0 

C valid  fur  EMACH  GTR  0.9 

6 EHg  = t MACH-().9 
HOLn=l  .-EXP  (-5.<>EM9) 

XKL  = n.550  7 + 0.  A<>hOLD 
XKL=XKL*6.6"SYaw 
XCP  = 0 .237* 1 .57«H0LD 

C VALlt'  FUR  LOACH  GTl>!  1.0 

7 XKS=-A.0  + 1 .7(3«  (EHaCH-1  . ) 

C valid  FOR  LOACH  OTR  1.1 

XKDYA'.v  = l .5  + 2. 3R»LXR  (-2. 720  (EMACH-1 
C VALiL'  FOR  LmaCH  GTR  1.3 

XKH=3.7 

C VALI''  FOR  ALl.  EnaCH 

9  SMAR(':  = XCP-XCG 
XKM  = -XKL<^SMARG 
IF (NaRTRL  .NE.O)  GO  TO  51 
RETURN 

10  XKA=n,nn5B 
XKDYaW=  1.5 

11  XKL  = 0.62-0.077<:-EHACH 
XkL=XKL+A.3oSYaW 
XCP=1 .2-1 .07«EHACH 

12  XKS=-A.O 

13  XM-|=n.71+2.3<tE*''ACH 
GO  TO  9 

20  EMB=FMaCH-O.B 

XKA=0.0038+0.nn?oEXP(-l .5"EM3) 
XKUYaw=1 .5+2.50SIN (6.263»EH8) 

GO  TO  11 

30  EMP.=EMaCH-O.B 

XKA=0.0030  + 0.nn?»EXP  (-1 .5«ENFi) 

XK['YaW=1 .5  + 2.5«SIiJ  (6.2H30EM8) 

CM9=EMaCH-0.9 

HOLn=l .-EXP (-5.*rM9) 

XKL=n.5507+0.A«HOL0 

XKL=XKL*5.5»SYAW 

XCP=o. 237+1 .57OH0LD 

GO  TO  12 

35  EM8=EMaCH-0,8 


O”! 10200 
00110300 
OOllOAOO 
00110500 
00110600 
00110700 
00110800 
00110900 
00111000 
00111100 
0011 1200 
0011 1300 
001 1 lAOO 

8))  00111500 

00111600 

00111700 

00111800 

00111900 

00112000 

00112100 

00112200 

00112300 

00112900 

1 ) ) 001  12500 

00112600 
00112700 
00112800 
00112900 
00113000 
00113100 
00113200 
00113300 
00113900 
00113500 
001 13600 
00113700 
00113800 
00113900 
00119000 
00119100 
00119200 
00119300 
00119900 
00119500 
00119600 
001 19700 
00119800 
00119900 
00115000 
00115100 
00115200 
00115300 
00115900 
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XKA=0.003B+0. on2«EXP(-l .b^EMS) 

XKDYAWrl ,5+2,5*SIN(6,2H3«EM8) 
EM9=EmaCH-0,9 
HOl.U=l  .-LAP  (-5.  *FM9) 
XKL?0.'5507+0.4<hiolO 
XKL=XKL*5.b*SYAA 
XCP=0.237+1 .57»h0lD 
XKS=-5.7B*1 .7B«EmaCH 
GO  TO  13 

40  EMS=FMaCH-O.H 

XKA=O.003rt*0.no2"EXP(-l 

XKOYAWrl  , 5*2. 36*EXP  (-2,72*  (E‘^ACH-1  , 1 ) ) 

LK9=E;'1aCH-0.9 

H0L0=l.-EXP(-5.*eM9) 

XKL=0.b507*0 .4*H0L0 
XKL  = XKL*6. 6*5YA»«. 

XCP=0.237*1  .r>7*HnLO 
XKS  = -5. 7tU  1 . 70*EMACH 
GO  TO  13 

45  XKA=0.0038*0.002*EXP(-1 ,5* (EMACH-0,«) ) 
GO  TO  h 
END 


00115500 

00115600 

00115700 

00115800 

00115900 

00116000 

00116100 

00116200 

00116300 

00116400 

00116500 

00116600 

00116700 

00116800 

00116900 

00117000 

00117100 

00117200 

00117300 

00117400 

00117500 

00117600 
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SUBROUTINE  RUNGEl (V,W»NEQiNORD»DIFEQ) 
DIMENSION  V(12)  ,W(4a) 

NV=NEQ*NORU 

N=MV4.NEU 

RETURN 

ENTRY  RUNGE2{T,DT) 

DT2=()T«.5 

DT6=DT/6. 

DO  1 1=1, N 

1 W{I)=V{1) 

DO  2 J=1 ,3 
NJM=MEQ  + N<>  { J-1 ) 

JDECK=N*J 

IF  (J-3)3,'+,4 
3 OTW=nT2 
GO  TO  5 
U DTV(=OT 
STvJ  = T^OTw 
DO  6 I=1,NV 
K=l+JOECK 
L=I^NJM 

W{K)=W{I)  ♦W{L)<>DTW 

6 V(I)=iV{K) 

call  DIFEO {Tw,V, J) 

DO  2 I=1,N 
K=I^JOECK 

2 W(K)=V(I) 

DO  7 1=1, NV 
Kl=l+NEO 
K2=K  1 4-0 
K3=K2+N 
KA=K3*N 

7 V{I)=W{l)+f)T6»{U{Kl)+2.*{W(K2)^W(K3)) 
T = TW 

CALL  DlFEQ {T,V,4) 

return 

End 


00117700 
00117800 
00117900 
00118000 
00118100 
00118200 
00118300 
00118400 
00118500 
00118600 
00118700 
00118600 
00118900 
00119000 
00119100 
00119200 
00119300 
00119400 
00119500 
00119600 
00119700 
00119800 
00119900 
00120000 
00120100 
00120200 
00120300 
00120400 
00120500 
00120600 
00120700 
00120600 
00120900 

W(K4))  00121000 

00121100 

00121200 

00121300 

00121400 


t 
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FUNCTION  VWC(Y) 

IF(Y.GE.6700.)GOTO3 

VWC=5.1816*4.972E-5*Y+1.3A9AE-7*Y*Y 

RETURN 

3 VWC=10.058*3.962A*COS(.42946E-3*(Y-9448,8) ) 

IFtY.GT. 13700.)  WRITEtGtl) 

RETURN 

1 FORMAT (28H  ALTITUDE  ABOVE  13700  METERS) 

END 


00121500 

00121600 

00121700 

00121800 

00121900 

00122000 

00122100 

00122200 

00122300 
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SUBROUT  I NE  BURN ( T IME » XMASS , THRUST ) 

subroutine  COMPUTES  PROJECTILE  MASS  IN  POUNDS  MASS  AND 

rocket  thrust  in  POUNDS  FORCE 

T0=TIME  at  InHICH  burning  commences 

EM0=1NITIAL  I-iASS.  lbm 

EMB=PURNT  MASS,  LBM 

TIME=TIME  after  launch  ,sec 
SPT=SPECIFIC  IMPULSE,  LBF/LBM/SEC 

fc=constani  nominal  thrust  level,  LBF 

It<URN=  INDICATOR  OF  COMMENCEMENT  OF  dURNI  NG  ( I DURN=  1 ) 
DELT  = RISE  time  OF  THRUST  — ASSUMED  EU'JAL  TO  DECAY  TIME  , 
BRATF  = FC/SPI=RIIRNING  RATE,  LBM/SEC 
TO= (EM0-EMB)/8RaTE=EFFECTIVE  BURNING  TIME,  SEC 


SEC 


COMMON  emo,emb,spi ,fc,brate,delt,to»tb, ISW, V, theta, FFCTR,CaLSO, 

1 VW,VCW,ALT ,R, iemd,cmach,rlynld,resib,cal,dlong, iopty,yaw,amom, 

2 BMOf!,PSI  ,'vTARFA,  ISEP,NABLE,  IGLIDE,  ING,TENABL,  AVTHD,CDEFL, 

3 MMCSW,TM0MMC 


00122A00 

00122500 

00122600 

00122700 

00122800 

00122900 

00123000 

00123100 

00123200 

00123300 

00123AO0 

00123500 

00123600 

00123700 

00123800 

00123900 

0012A000 

0012A100 


TO)  GO  TO  1 

0012A200 

TO+DELT)  GO 

TO  2 

0012A300 

0012AA00 

T2)  GO  TO  3 

0012A500 

T2  + l'LLT)  GO 

TO  4 

0012A600 

T2=T0*TB 
IF  (TIME .LE 
IF (TIME.LT 

THf!UST  = 0.  0012A700 

XMA5S=EMB  0012A800 

RETURN  0012A900 

XmaSS=EMO  00125000 

TIIRU5T  = 0.  00125100 

RETURfJ  00125200 

XMASS=EMO-RRaTFo  (TlME-TO)*<>2/(2.«nELT ) 00125300 

THRUST = ( T IME-TO)/DELT*FC  00125A00 

RETURN  00125500 

XMA5S  = EM0-r:RATE«  (TIME-TO-DELT/2.  ) 00125600 

THRUST=FC  00125700 

RETURN  00125800 

XMA5S  = EM0-BRATF» ( (TB-DELT/2, ) ♦ (TIME-T2) * ( 1 . - ( T IME-T2) /DELT/2. ) ) 00125900 

THRUST=FCo (1 .- ( T IME-T2) /DELT)  00126000 

RETURN  00126100 

END  00126200 
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FUNCTION  I8URN(TIME»ALT»QE»VEL0) 

FUNCTION  PRODUCES  INDICATION  OF  COMMENCEMENT  OF 
BURNING — I BURNS  1. 

IBUKNiO  UNTIL  BURN  BEGINS 

USER  CHOOSES  TO  FROM  CURREtjT  TIME  OR  ALT  (ALTITUDE)  OR 
QE  (LOCAL  QUADRANT  ELEVATION)  OR  VELO  (VELOCITY) 

COMMON  /SRCOM/TM 
DATA  ALTMAX/OOOOn./tVMIN/O.O/ 

IF (7IMF.6E.TM)  GO  TO  3 
IF (ALT.GE.ALTMAX)  go  to  3 
IF (VELO.LE. VMIN)  GO  TO  3 
C IF (OE.GT.45.0)  go  TO  2 

IBURr.'=0 
return 

C 2 IF (OE.LT.AG.n)  60  TO  3 

C A IBUR!I  = 0 

C RETURN 

3 IUURn=1 

RETURN 
LNO 


00126300 

00126A00 

00126500 

00126600 

00126700 

00126800 

00126900 

00127000 

00127100 

00127200 

00127300 

00127A00 

00127500 

00127600 

00127700 

00127800 

00127900 

00128000 

00128100 

00128200 

00128300 
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BLOCK  data  001?8AOO 

DIMENSION  TK/i  ( 1 1 ) ,TKDYAW  ( 1 1 ) »TKL  ( II  > *TKM(  1 1 ) *TKF  (ID*  00128500 

1 TKT(ll) *TKM(11) *TK5(11) *TCP(11) *TMACh(11)  00128600 

DIMENSION  TF^'ACH(4)  ,ALTRmm(A)  *AALTRM(4)  ,TCN(A*4)  *V«TCN(A)  00128700 

DIMENSION  TCaDl(4),TCA02(A) , IDEL (A)*7aT(4,A)*WTAT(A)  00128800 

COMMnN/COFCOr’/XCG*SMAKG*EM*  THIL)*PRnU  , ALTR  IM , CNATRM * GL  IDE  , EPSTHE  00128900 

1 ,5TaFAC*yAwNU, TKA,TkOYAa  * TKL*TKM*TKF ,TKT*  TKH*TKS*TCP*TMACH*  00129000 

2 MARTtH,  *TFMACH,  ALTRMM*AALTKM,TCf4*WTCN,TCADl  , T CA  D2  * TDEL  * T AT  * W T AT  00129100 

CO***  OLFADLT  VALUES  FOR  AFRu  UAlA  USED  IN  GLIDE  CALCULATIONS  00129200 

data  TFMACH( 1 ) /0.50/*  TFMACH (2) /0.8  0/*  TFMACH (3) /0.90/*TFMACH (4 ) / 1 . /O 0 1 2930 0 
1*ALTPMM(1)/14.30/*ALTRMM(2)/14.00/*ALTRMM(3)/13.B0/*ALTRMM(4)/13,60 0129400 
2 0/* AALTRM(I) /0.0/*AALTRM(2) /5.0/*AALTRM(3) /in.0/*AALTRM(4)  00129500 

3/15.0/,TCN(1,1)/0.0/*TCN(1,2)/1.10/*TcN(1,3)/2.10/*TCN(1*4)/3.10/*00129600 
4TCN (2, 1 ) /0.0/*TCN (2*2) /I .20/*TCN(2,3) /2. 10/*TCN(2*4) /2.90/,  00129700 

5TCN(3* 1) /0.0/*TCN(3,2) /I .30/*TCN(3,3) /2 . 20/ * TCN ( 3 , 4 ) /3 . 00/ , 00129800 

6TCN(4, l)/0.0/*TCN(4,2) /I . 30 / * T CN ( 4 , 3 ) /2 . 50/ * TCN ( 4 , 4 ) /3 . 7 0/  00129900 

DATA  TCAD 1/4, 0*5. U*5.75*fc>.50/*TCAD2/b, 0*5. 0*4. 5*4.0/*  00130000 

1 TOEL/n. *5.0*ln.0*15.0/*TAT/0.*0.*0**0.*7.4*7.2*7.4,7.1*  00130100 

2 12.1*11.8*11.6*11.3*17.7*17.2*17.1*17.1/  00130200 

end  00130300 
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MACH  NO  COEF  D«AG  DkaG  IMCR 
0.0  0.3250  0.0350 
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DRSAR-SAM 


I * JUN  1977 

MEMORANDUM  FOR  RECORD 

SUBJECT:  Information  for  SHAPE  Technical  Center  Relative  to  Z0T.14/ 
Z0T.15 


1.  Reference  is  made  to  Technical  Note,  DRSAR/SA/N-58,  Description  of 
a Computer  Program  (Z0T.14)  for  Guidance  Simulation  of  Cannon-Launched 
Guided  Projectiles  (AD  #A036663),  Jan  77. 

2.  DRSAR-SA  is  preparing,  for  transmittal  to  SHAPE  Technical  Center,  a 
set  of  data  including  the  source  program  for  Z0T.15,  the  guidance  simu- 
lation program  currently  used  by  DRSAR-SAM.  This  memorandum  provides  a 
guide  to  use  of  this  program. 

3.  User  information  required  for  use  of  this  program  may  be  classified 
as  (a)  information  relative  to  conversion  of  the  source  from  IBM  S-360 
extended  FORTRAN  to  ANS  FORTRAN  for  UNIVAC,  and  (b)  methodological  back- 
ground with  input/output  descriptions. 

4.  Relative  to  the  first  category  (conversion),  the  following  items  need 
be  pointed  out: 

a.  Precision  need  not  be  double  in  the  UNIVAC  version  due  to  longer 
word  length.  IMPLICIT  statements  are  to  be  deleted,  and  variables  declared 
simply  REAL  or  INTEGER. 

b.  Certain  I/O  formats  will  require  change — fields  specified  as  D 
should  be  changed  to  E,  literal  (Hollerith)  fields  should  be  rewritten  in 
H-format  (the  "quote"  format  has  been  used),  and  the  lengths  of  some  A- 
formatted  fields  may  require  change  (A8  has  been  used). 

c.  Numeric  literals  (constants)  written  D- format  throughout  the  pro- 
gram should  be  changed  to  E-  or  F-format  as  required. 

d.  Multiple-entry  subroutines  may  not  be  allowed.  If  not,  such  a 
subroutine  may  be  separated  into  several  subroutines  sharing  common  stor- 
age, or  alternatively,  may  be  converted  to  a single-entry  subroutine  with 
a parameter  in  the  call  list  specifying  which  segment  is  to  be  performed 
(DODGER  is  an  example). 

e.  Subroutine  RANDMM  is  the  pseudo-random  number  generator  used  in 
this  program.  It  is  a slight  re-code  of  RANDU  of  the  IBM  S-360  Scientific 
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» 


f 2 JUN  1977 


DRSAR-SAM 


SUBJECT:  Information  for  SHAPE  Technical  Center  Relative  to  Z0T.14/ 
Z0T.15 

Subroutine  Package  (SSP).  This  code  is  applicable  only  to  the  360  and 
must  be  replaced  by  code  applicable  to  whatever  machine  is  to  be  used, 

5.  Relative  to  the  second  category  (methodology  and  I/O),  the  referenced 
technical  note  provides  the  basis  for  understanding  the  present  program. 
The  key  difference  between  the  two  programs  is  in  subroutine  TRACK,  which 
was  totally  changed,  in  order  to  provide  a more  faithful  representation 
of  the  logical  states  of  the  seeker.  Other  changes  are  generally  of  the 
nature  of  an  added  option.  These  are  discussed  in  the  input  guide, 
attached  (Incl  1),  which  also  describes  the  input  data  required  to  run 
the  program. 


1 Incl 


RICHARD  D.  HEIDER 


as 


Operations  Research  Analyst 
Methodology  Division 
Systems  Analysis  Directorate 
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USER'S  INPUT  GUIDE  TO  ZOT.15 


The  following  input  list  is  provided  as  an  amendment  to  the  correspond- 
ing list  found  in  Appendix  A of  the  tech  note  referenced  in  the  body  of 
this  memorandum. 

CARD  1:  Identical  to  CARD  1 of  Z0T.14 

CARD  2:  Identical  to  CARD  2 of  Z0T.14 

CARD  3:  Identical  to  CARD  3 of  Z0T.14 

CARD  4,  required,  format  (8F10.0),  contents  by  field: 

(1)  THETGY  elevation  angle  of  gyro  (optical)  axis  [deg],  ignored 

unless  IGLIDE  is  equal  to  2 (see  CARD  9) 

Remaining  fields  are  not  used. 

CARD  5:  First  seven  fields  identical  to  corresponding  fields  of  CARD  4 
of  Z0T.14.  Field  (8)  contains 

(8)  CANT  gyro  (optical)  axis  cant  or  lookdown  angle  [deg],  appli- 

cable to  ballistic  mode.  (Do  not  use  this  variable  to 
account  for  hangoff  as  suggested  with  Z0T.14,  since 
this  effect  is  now  properly  modeled.) 

CARD  6;  Identical  to  CARD  5 of  Z0T.14 

CARD  7:  Identical  to  CARD  6 of  Z0T.14 


CARD  8,  required,  format  (8F10.0),  contents  by  field: 


(1) 

CSMAl 

as  in  Z0T.14 

(2) 

CSMA2 

as  in  Z0T.14 

(3) 

ROLRAT 

as  in  Z0T.14 

(4) 

CTERM 

control  parameter— specify  0.0  to  terminate  the  run 
upon  any  impact  for  which  acquisition  had  failed  to 
occur;  specify  1.0  to  continue  the  run  to  the  next 
replication  regardless. 

• 

(5) 

TMAX 

control  variab1e--specify  0.0  to  exercise  the  internally- 
defined  time  limit  (standard  with  Z0T.14);  specify  a posi' 
tive  clock  time  [sec]  (greater  than  maximum  expected  im- 
pact time)  to  exercise  option  of  user-specified  time  limi 

(6) 

KLDOT 

as  in  Z0T.14  except  that  it  is  now  applicable  to  both 
glide  and  ballistic  modes. 

1 
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(7)  GSBAR 

range  of  gimbal  angle  (abs.  val.)  within  which  the 
gimbal  saver  output  is  zero  [deg] 

(8)  KGS 

gain  of  gimbal  saver  [dimensionless] 

CARD  9,  required,  format  (7F10.0,  212,  213),  contents  by  field: 


(1-7) 

identical  to  corresponding  fields  of  CARD  8 of  Z0T,14 

(8)  IGLIDE 

projectile  mode  control  parameter--0  for  ballistic,  1 
for  glide  trim,  2 for  attitude-hold  with  user-specified 
attitude.  Note:  1 and  2 are  both  glide  options— 

IGLIDE  = 1 causes  projectile  and  gyro  attitudes  to  be 
internally  defined  such  that  the  initial  acceleration 
of  the  projectile  lateral  to  the  velocity  is  zero 
(standard  with  Z0T.14),  while  IGLIDE  = 2 causes  pro- 
jectile and  gyro  attitude  to  be  as  specified  by  THETGY 
on  CARD  4 

(9)  ICAGE 

gyro  cage  control  parameter  for  action  of  gyro  after 
loss  of  acquisition-- ICAGE  = 0 for  free  gyro  or  1 for 
caged  gyro 

(10)  NLAST 

output  option--NLAST  = 0 to  ignore,  or  1 to  100  to  aver- 
age the  positions  of  the  last  NLAST  spots  and  compute 
and  report  the  miss  distance  with  respect  to  this  aver- 
age position  in  addition  to  the  miss  distance  with 
respect  to  the  nominal  aimpoint  (pitch,  yaw,  and  modulus) 

(11)  NEWAER 

new  aero  option— specify  0 to  use  the  aero  tables  cur- 
rently stored  (defined  by  BLOCK  DATA,  if  this  is  the  first 
run  of  the  job,  or  by  input  from  a previous  run),  or  spec- 
ify 1 if  any  changes  are  to  be  made. 

CARD  10,  optional,  format  (1615),  contents  by  field: 

(1)  N1  option  parameter  for  vector  VMACH--N1  = 0 for  no  change. 


N1  = 1 thru  7 for  a changed  vector  VMACH  having  NMACH  = 

N1  elements 

(2)  N2 

similar  parameter  for  VMCAO,  NMCAO 

(3)  N3 

similar  parameter  for  VDELTA,  NDELTA 

(4)  N4 

similar  parameter  for  VALFHA,  NALPHA 

(5)  11 

option  parameter  for  table  TKA — 11  = 0 for  no  change, 

1 for  change 
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(6) 

12 

similar  parameter  for  TCNA 

(7) 

13 

similar  parameter  for  TCND 

(8) 

14 

similar  parameter  for  TCAD 

(9) 

15 

similar  parameter  for  TXSM 

(10) 

16 

similar  parameter  for  TCAO 

Remaining  fields  are  not  used. 

CARDS  11-i,  optional,  format  (8F10.0),  contents:  Each  vector  or  table  for 
which  a change  is  required  is  punched  in  the  order  of  mention  in  the 
above  description  of  CARD  10;  each  vector  or  table  begins  a new  card;  for 
tables  TKA  and  TCNA,  rnach  number  varies  before  varying  alpha  or  delta,  as 
appropriate.  See  the  listing  of  subroutine  AEROIN  for  further  information 
on  using  this  option. 

CARD  12:  Identical  to  CARD  9 of  Z0T.14 

CARD  13:  Identical  to  CARD  10  of  Z0T.14 

CARD  14:  Identical  to  CARD  11  or  Z0T.14 


CARD 

15,  required 

, format  (8F10.0),  contents  by  field: 

(1) 

BMDIVG 

as  in  Z0T.14 

(2) 

AGCLD 

same  as  DYRANG  of  Z0T.14 

(3-7) 

identical  to  corresponding  fields  of  CARD  12  of  Z0T.14 

(8) 

KC 

gyro  electrical  cage  gain  [sec”^] 

CARD 

16,  required 

, format  (8F10.0),  contents  by  field: 

(1) 

AGCSR 

max  rate  [db/pulse]  at  which  detector  automatic  gain 
control  (AGC)  can  shift  (up  or  down) 

(2) 

WLE 

nominal  width  of  loading  edge  of  detector  glimpse  gate 
[microsec].  This  is  the  time  between  opening  of  the 
gate  and  expected  arrival  of  the  pulse  from  the  target. 

(3) 

WGATE(l) 

width  of  gate  during  correlation  sequence  [microsec]. 
Time  from  opening  to  closing  of  gate  if  no  pulse  is 
received. 

(4) 

WGATE(2) 

width  of  gate  after  correlation  is  established  (during 
tracking). 
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(5) 

WTRUNC 

time  from  reception  of  first  pulse  in  gate  to  gate 
closure  [microsec].  Supersedes  closure  time  scheduled 
by  WGATE 

(6) 

GATESR 

max  rate  [microsec/pulse  interval]  at  which  the  leading 
edge  (opening)  of  the  glimpse  gate  can  be  shifted  rela- 
tive to  nominal  pulse  period. 

(7) 

SLOPE 

slope  of  intrapulse  dynamic  threshold  [db/microsec] 

(8) 

not  used 

CARD 

17,  required 

, format  (8F10.0),  contents  by  field: 

(1-5) 

(6)  XINCR 


(7)  ZINCR 


(8)  CRELOC 
CARD  18:  Identical  to  CARD  13-2  of  Z0T.14 
CARD  19:  Identical  to  CARD  14  of  Z0T.14 
CARDS  20:  Identical  to  CARDS  15  of  Z0T.14 


CARD  21 
field: 

, optional 

, format  (8F10.0)— used  only  when  MODESM  = 3— contents  by 

(1) 

RDB 

distance  [m]  from  designator  to  nearest  point  of  back- 
ground, which  is  modeled  as  a vertical  plane  of  infinite 

extent 

(2) 

AZDB 

azimuth  [deg]  from  designator  to  the  same  point 

(3) 

AZDES 

azimuthal  direction  [deg]  of  designator's  beam 

(4) 

ELDES 

elevation  [deg]  of  the  same 

(5) 

REFL2 

reflectivity  of  background 

identical  to  corresponding  fields  of  CARD  13-1  of  Z0T.14 

if  CRELOC  = 0.0,  XINCR  is  the  x-wise  increment  as  de- 
fined for  Z0T.14;  if  CRELOC  = 1.0,  XINCR  is  the  standard 
deviation  of  a circular-normal  distribution  of  target 
locations  [m;ft] 

if  CRELOC  = 0.0,  ZINCR  is  as  defined  for  Z0T.14;  if 
CRELOC  = 1.0,  XINCR  is  a positive  odd  integer  value  used 
as  the  seed  for  random  target  relocation  (number  of 
digits  should  be  less  than  the  number  of  significant 
decimal  places  allowed  by  machine  word  size  for  a real 
variable  of  the  precision  used). 

option  code  as  explained  above.  , 


70 


CARD  22:  Identical  to  CARD  16  of  Z0T.14 


mcdatP*  ^'ormat  (9X,5F8.2)~used  only  if  MODESM 

NSPOTS  cards  required— contents  by  field: 


2 or  3,  then 


(1) 

YV(1) 

as  in  Z0T.14 

(2) 

YV(2) 

as  in  Z0T.14 

(3) 

YV(3) 

as  in  Z0T.14 

(4) 

CROSS 

same  as  BRITE  of  Z0T.14 

(5) 

PCTEBG 

percent  of  beam  energy  spilling  over  onto  background. 

This  completes  the  data  deck. 


Next  page  Is  blaiik. 
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MEMORANDUM  FOR  RECORD 

SUBJECT;  Computer  Simulation  Study  of  the  Relationship  of  the  COPPERHEAD 
Footprint  to  Ceiling  and  Gun-to-Target  Range 


1 , References : 

a.  Conversation  between  Mr.  George  Schlenker,  DRSAR-SAM,  and 
COL  Robert  Nulk,  DRCPM-CAWS-GP,  24  Mar  77,  subject  as  above. 

b.  Memorandum  for  Record,  DRSAR-SAM,  31  Jul  76,  subject:  Computer 
Simulation  Study  of  COPPERHEAD  (CLGP)  for  Guidance  Accuracy  and  Footprint. 

c.  Technical  Report,  Rodman  Lab,  R-TR-77-007,  Jan  77,  title:  a' 
Comprehensive  Digital  Flight  Simulation  of  the  Cannon  Launched  Guided 
Projectile. 

2,  Introduction. 

This  memorandum  documents  a study  performed  by  the  undersigned,  with 
the  guidance  of  G.  Schlenker,  during  Apr-Jun  77  at  the  request  of  COL  Nulk 
(Ref  la),  to  further  explore  the  dependence  of  footprint  size  upon  the 
ceiling  and  th6'’gun-to-target  range  (GTR).  The  most  recent  previously- 
developed  estimates  of  footprints  for  COPPERHEAD  (Ref  lb)  included  two 
levels  of  ceiling  (3000  and  2000  feet)  and  three  GTR's  (6,  12,  18  kilo- 
meters). The  intent  of  this  study  was  to  generate  footprints  for 
unlimited  ceiling  and  for  ceilings  from  3000  feet  to  zero,  by  decrements 
of  500  feet,  for  the  same  set  of  three  GTR's.  At  the  same  time,  certain 
of  the  projectile's  design  parameters  were  to  be  updated,  in  order  that 
the  estimates  might  be  as  current  as  possible. 

3,  Parameter  Changes. 

Many  of  the  seeker  and  autopilot  parameters  were  altered  since  the 
previous  study  (Ref  lb).  Those  whic_h_one  might  expect  to  affect  the  foot- 
print include:  Kg,  K^,  K^,  Kgg  and  GS.*  In  addition,  flight  tests  per- 
formed in  Mar  77  indicated  a larger  drag  coefficient  than  previously 
estimated. 


*These  parameters  are; 

Kg  gain  of  g-bias  computation  circuit  K^g  gain  of  gimbal  saver 

Kg' gain  of  gyro  electrical  cage  circuit  5^  threshold  of  gimbal  saver 

K^  gain  of  attitude-hold  circuit 
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4.  Procedures. 

a.  Zoning  Solutions 

The  EXBAL  exterior  ballistic  program  was  used  to  generate  unguided 
trajectories  for  the  projectile  using  new  aerodynamic  drag  estimates. 

From  this  set,  appropriate  trajectories  for  the  three  GTR's  were  selected 
(Table  1). 

b.  Footprint  Definition 

The  Z0T.15  guidance  simulation*  was  then  used  to  generate  the  foot- 
prints shown  in  Figures  1-3  and  listed  in  Table  2. 

c.  Experiments  for  Stretch  Range  Dependence 

Further  Z0T.15  experiments  were  performed  to  define  the  stretch  range 
of  COPPERHEAD,  given  the  stated  12-km  trajectory,  as  a function  of  ceiling 
and  meteorological  visibility.  Results  are  displayed  graphically  in 
Figures  4 and  5. 

Figure  4 shows  the  dependence  upon  visibility,  and  Figure  5 shows  the 
dependence  upon  ceiling.  For  any  given  combination  of  ceiling  and  visi- 
bility, the  stretch  range  is  the  lesser  of  the  two  stretch  ranges  indi- 
cated from  these  figures. 

5.  Analysis  of  Results. 

Generally,  the  dependence  of  the  footprint  upon  ceiling  and  upon  GTR 
is  similar  to  that  seen  previously  (Ref  lb).  For  low  ceilings,  however,  the 
size  of  the  footprint  shrinks  extremely  rapidly,  due  to  the  very  short  time 
available  for  proportional-navigation  (PN)  guidance.  Note  that  the  foot- 
print has,  for  all  practical  purposes,  vanished  at  the  1000-foot  level.  In 
fact,  before  reaching  the  500-foot  level,  the  footprint  vanishes  totally, 
because  the  time  remaining  to  impact  is  insufficient  for  the  seeker  to 
sequence  to  PN  guidance  or  to  arm  the  warhead.  The  fly-to-seeker  guidance 
(FTS)  is  inadequate  for  purposes  of  hitting  the  target;  an  interval  of 
several  seconds  of  PN  guidance  is  essential. 

The  footprints  generated  during  these  experiments  agreed  well  with 
those  of  Ref  lb  except  in-the  stretch,  which  is  significantly  reduced  in 


♦This  program  is  the  successor  to  Z0T.14,  documented  in  Technical  Note 
DRSAR/SA/N-58  (AD#  AO  36663),  January  1977,  Description  of  a Computer  Program 
(Z0T.14)  for  Guidance  Simulation  of  Cannon-Launched  Guided  Projectiles.  The 

supplementary  documentation  for  the  current  Z0T.15  is  contained  in  MFft, 

DRSAR-SAM,  22  Jun  77,  subject:  Information  for  SHAPE  Technical  Center  Rela- 
tive to  Z0T.14/Z0T.15. 
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parh  case  Further  experiments  were  performed  using  Z0T.15  to  identify  the 

reduce^  the  navigation  gain  whenever  the  projectile  axis  is  steerea  mo 
than  GS  away  from  the  gyro  axis. 

6.  Verification  of  Gimbal -Saver  Effect. 

It  was  desired  to  cross-check  between  simulation  models  9'^ 
rtf  the  effect  of  the  gimbal -saver  upon  the  stretch  range  as  indicate  y 
5nT  1^  Therefor^  a series  of  runs  of  the  Rodman  six-degree-of-freedom  _ 
mSde^(Ref^^c)  were  made  for  DRSAR-SAM  for  a siightiy  different  12-km  nominal 
GTR  traiectory  (Rodman  aero  being  slightly  different  from  that  used  in  . 
p^L^ding  portion  of  this  study).  Z0T.15  matching  runs  7’" 
exact-match  inputs,  with  the  results  displayeo  in  terms  of  miss  distance  vs 
range  displayed  in  Figure  7.  These  results  do  verity  the  existence 
Sfea  of  the  Qimbal-sever,  but  there  is  a disagreement  between  the  two  ^ 
simulation  models  as  to  the  stretch  point,  which  is  presently  unexplaineo. 
However,  there  is  adequate  agreement  between  these  models  as  to  the  tuck 
point.  Further  analysis*  of  differences  is  requireo. 

7.  Conclusions  and  Recoirmiendations . 

a For  ceiVifigs  much  below  2000  feet,  an  adequate  footprint  is  not 
achieOable  u"ng  thi  desired  20-degree  glide;  8 ide  might  p- 

vide  the  reauired  footprint,  but  one  must  accept  decreaseo  lethali  y 
that  apprSrchTs  taken^  TrJde-off  studies  along  this  Ime  are  reco-enoed. 


b. 


The  effect  of  the  gimbal-saver  requires  further  study,  preferably 
using  all  Laluble  flight  models.  Eased  upon  the  results  of  this  study, 
I recommend  a review  of  the  choice  of  the  parameter  values  of  the  9^7!^ 
Lvfr  wUh  a viL  to  restoring  the  gimbal-saver  to  its  original  design  or 

its  possible  elimination. 


*Additional  runs  would  be  required  of  the  Rodm.an  tJp^transfer- 

turnaround  would  be  excessive  due  to  the  demands  im, posed  by  the  trans  e 

of-function  process  currently  underway. 


9 Incl 

1-7.  Figure's 

8.  Table  1 

9.  Table  2 


RICHARD  HEIDER 
Operations  Research  Analyst 
Methodology  Division_ 
Systemis  Ainalysis  Division 
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Figure  7.  Comparison  of  Z0T.15  and  Rodman  6-DOF  Results  for  Ginibal  Saver  Check  Runs 
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QUADRANT  ELEVATION  26,0.  26.75.  4A, 

GLIDE  ANGLE  (nominal)  (bal,).  -20.  -20 

PROJECTILE  DESIGN  - CURRENT  (APR  77)  EXCEPT  DRAG  INCREASED  OVER 
WIND  TUNNEL  RESULTS  PER  MAR  77  TEST  FIRINGS 
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® 7 m 1977 


MEMORANDUM  FOR  RECORD 

SUBJECT:  Scranton  Billet  Crane  Simulation 


1.  After  the  Scranton  billet  crane  simulation  project  was  completed  in 
August  1976,  Chamberlain  Manufacturing  Corporation  obtained  a copy  of  the 
report  and  wrote  a letter  to  ARRCOM  listing  their  objections  to  the  assump- 
tions made  in  the  study.  The  letter  was  forwarded  to  this  office  for  comments. 
Several  of  their  objections  were  judged  valid.  At  the  request  of  DRSAR-IMB 
changes  were  made  in  the  program  and  the  simulation  was  rerun  to  answer  the 
valid  objections.  This  memo  details  the  changes  made  and  the  analysis  per- 
formed. 

2.  In  summary,  the  simulation  results  indicate  that  the  200  feet-per^ minute 
maximum  velocity  recommended  by  the  Corps  of  Engineers  for  the  Scranton  AAP 
billet  crane  is  adequate. 


3.  The  following  new  information  was  supplied  by  DRSAR-IMB: 


billet 

number 

si  ze 

mult  wt. 

mults/mo. 

5 1/4" 

107  lbs 

147,000 

6 3/4" 

172  lbs 

21 ,000 

7 3/8" 

220  lbs 

63,000 

b.  An  automatic  squaring  table  is  being  procured. 

c.  The  crane  has  a hoist  speed  of  90  ft/mi n. 

d.  The  maximum  stacking  height  is  15  ft. 

4.  The  following  information  was  obtained  from  Chamberlain's  letter  or 
from  phone  conversations  with  Mr.  Bernie  White,  Scranton  AAP: 


billet 

size 

billet 

weight 

number 

billets/ht 

no.  billets 
per  charge 

5 1/4" 

1874  lbs 

170 

10 

6 3/4" 

3091  lbs 

103 

8 

7 3/8" 

3698  lbs 

87 

6 
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b.  The  average  load  per  rail  car  Is  60  tons. 

c.  There  are  500  work  hours  per  month.  Including  coffee  and  shift  breaks. 

d.  Rail  car  delivery  Is  restricted  to  16  hours  per  day. 

e.  Breaker  line  capacity  Is  340  mults/llne/hr. 

5.  The  following  Is  an  explanation  of  the  changes  made  In  the  assumptions 
and  data  Inputs  used  by  the  simulation.  Most  of  the  data  used  In  the  sim- 
ulation was  collected  at  the  time  of  the  original  study. 

a.  Assumption  5.  Billets  are  stored  and  charged  into  the  feeders  In 
heats  which  are  assumed  to  be  groups  of  170,  103,  and  87  each  for  5 1/4".  • 

6 3/4",  and  7 3/8"  billets,  respectively. 

Assumption  12.  The  crane  operator  (las  sufficient  skill  and  is 
allowed  to  operate  the  crane  at  maximum  speed.  He  could,  thus,  begin  x 
and  y movement  of  the  crane  simultaneously,  after  lifting  the  load  above 
the  rail  car  sides  or  bay  stock-piles. 

c.  The  random  number  generator  has  been 'changed  to  use  a separate 
string  of  random  numbers  for  each  stochastic  item.  Therefore,  the  first 
six  Input  cards  contain  42  seeds  In  an  8110  format. 

d.  I tern  1 , Chamberlain  states  that  under  mobilization  conditions 
they  will  work  500  hours  per  month,  including  coffee  and  shift  breaks. 

Since  there  are  45  minutes  of  breaks  plus  35  minutes  for  lunch  per  8 hour 
shift,  the  actual  number  of  work  hours  per  month  is: 


The  mean  time  between  charges  for  the  5 1/4"  billets  is  determined  by: 


number  of  mults  per  billet 


1874  Ibs/billet 

~1V/  Ibs/nwlt 

17  mults/blllets 


mean  time  between  charges  ■ 

(17  mults/bllletino  billets 


* 31.19  min/chg 


86 


DRSAR-SAL 

SUBJECT:  Scranton  Billet  Crane  Simulation 


7 JUN  W7 


To  create  some  manufacturing  variability,  It  will  be  assumed  that  this  31.19 
minutes  between  charges  Is  the  mean  of  a normal  distribution  having  95X  of 
Its  area  within  15X  of  Its  mean.  Thus,  the  standard  deviation  Is; 

- 2.39 

Cut-offs  are  set  at  3 std  dev  or  24.02  minutes  and  38.36  minutes. 

e.  Item  2.  The  distribution  of  the  time  between  charges  for  the  6 3/4" 
billet  table  Is  as  follows: 

number  of  mults  per  billet  * ^172"! b's^^l t^'^" 

* 18  mults/blllet 


mean  time  between  charges  = 

(18  mults/b111et)(8  b111ets/chq)(449.44  hr/mo)(60  mln/hr) 

(21 ,00Q  mults/mo) 

= 184.91  mln/chg 
_ _ (.15)(184.91) 

= 14.15 

limits  = 142.46  minutes,  277.36  minutes 

f.  Item  3.  The  distribution  of  the  time  between  charges  for  the  7 3/8" 
billet  table  Is  as  follows: 

number  of  mults  per  billet  - 

* 16  mults/blllet 


mean  time  between  charges  = 

(16  mults/blllet) (6  b111ets/chq)(449.44  hr/mo)(60  mln/hr) 

(63,000  mults/mo) 


■ 41.09  mln/chg 
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(■15j(41.09) 


3.14 


limits  = 31.67  minutes,  50.51  minutes 

g.  Item  4.  Since  the  arrival  of  heats  Is  not  affected  by  breaks.  It 
must  be  based  on  the  number  of  hours  per  month  that  the  plant  Is  operating: 

(5W  luncK)  ’ ”9-33  hr/mo  ' 

The  mean  time  between  arrivals  of  5 1/4"  heats  Is  determined  by: 


147,000  mults/mo 

17  mults/billets 


8647.06  billets/mo 


= 60.865  hts/HB 


(539.33  hrs/mo)(60  min/hr) 

(50.865  nts/mo)  ^ 


636.19  min/ht 


Since  most  arrivals  are  Poisson,  it  will  be  assumed  that  this  arrival  is 
also  Poisson  distributed.  Further,  since  these  arrivals  are  not  completely 
randotn,  cut-offs  at  3 standard  deviations  will  be  employed.  For  a Poisson 
distribution,  the  standard  deviation  is  equal  to  the  square  root  of  the 
mean.  Thus  the  limits  are  560.52  and  711.86.  Chamberlain  indicated  that, 
rail  cars  are  delivered  during  a 16  hour  period  each  day  and  heats  arriving 
during  the  remaining  8 hours  are  delayed.  The  program  has  been  changed  to 
delay  any  arrivals  which  occur  in  the  last  8 hours  of  any  24  hour  period. 


h.  Item  5.  Arrival  of  6 3/4"  heats. 


(21,000  mults/mo) 

(rS  muTt^TTTFtt 

I 


1166.67  billets/mo 


1166.67  mults/mo)  _ ,,  u*./-- 

T03  Vm^feTht) 
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(539.33  - 2856.91  mln/ht 

2856.91  +3(2856.91)^  - 2692.56,  3017.26 
Thus  the  limits  are  2692.96  and  3017.26. 

1.  Item  6.  Arrival  of  7 3/8"  heats. 

• = 3937.5  billets/™ 

(539.33  - 715  n.1„/ht 

715.00  +3(715)^  = 634.78,  795.22 
Thus  the  limits  are  634.78  and  795.22. 

j.  Item  11.  The  time  required  to  pick  billets  out  of  a rail  car  is 
entered  as  a triangular  distribution  having  a minimum  time  of  .40  minutes, 
a maximum  time  of  1.65  minutes,  and  a most-likely  time  of  .90  minutes. 

This  is  an  increase  of  9 seconds  to  insure  that  the  load  is  lifted  clear 
of  the  car  sides,  i.e.,  tolift  an  additional  13.5  feet. 

k.  Item  15.  The  number  of  5 1/4"  billets  picked  out  of  a rail  car 
per  unit  pick  is  entered  as  a triangular  distribution  having  1 as  the 
minimum,  12  as  the  maximum,  and  10  as  the  most-likely  number  of  billets. 

l.  Item  16.  The  number  of  6 3/4"  billets  picked  out  of  a rail  car 
is  entered  as  a triangular  distribution  with  1,  8,  and  6 as  the  minimum, 
maximum,  and  most-likely  number  of  billets,  respectively. 

m.  Item  17.  The  number  of  7 3/8"  billets  picked  out  of  a rail  car 
is  entered  as  a triangular  distribution  with  1,  8,  and  6 as  the  minimum, 
maximum,  and  most-likely  number  of  billets,  respectively. 

n.  Item  22.  The  time  required  to  pick  billets  off  of  the  storage 
pile  in  the  bays  was  entered  as  a triangular  distribution,  having  a minimum 
time  of  .40  minutes,  a maximum  time  of  1.15  minutes  and  a most-likely  time 
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cLfe  TirtercleIl:'of'?hrblrsl«L?'  ’ ’«« 

. °*  Items  23-25.  The  number  of  5 1/4".  6 3/4"  ^nH  7 

170  binllr^  5 1/4”  heat  Is  entered  es  . consUnt 

103  JilliH^  « 3/4-  heat  Is  entered  as  a constant 

87  bfilesf^  ^ 3/8"  heat  Is  entered  as  a constant 

39.  The  typical  number  of  5 1/4"  billets  on  a rail  ic 
entered  as  a constant  60  billets.  "lets  on  a rail  car  Is 

ente?;d  IrT^stl!!?  ® ^ '» 

entered  instil:?  If MUeEs!*^'’  ^ 

freas  5s1d°''Il^t!trel;aM IlSef “>  work 

The  simulation  was  ^un  th?ee  ?1mL  wUh  dliJ^LTr  ^ other  operations, 

generators.  Each  run  was  for  400  non  random  number 

and  billet  handling  Simis  f?om  ?Se  ?hrel  n!n  *1.,^?^  T'^f^^on  In  travel.  Idle, 
the  maximum  variation  was  less  than  dete^lned.  and  In  all  cases 

400.000  minutes  Is  iSnq  enoSah  Ps^ahiJch  indicates  that 

time  percentages.  The  result*;  nf  thlcoV^^  a steady  state  condition  for  the 
be  idle  at  least  332!  nf  tho  •n  runs  Indicate  that  the  crane  wouTd 

IncludI  cof||efsh1?t  an5  lunch  This  does  wt 

16.66%.  In  addUlJ;  ’theJe  ?s  a ol'n"  additional 

required  to  assure  that  thn  fPoH  P’^obablllty  that  the  maximum  lead  times 
8.92  minute!  and  8.8rmlSStL  ?o^1JrRS"/an®"r,^^^  9*97  minutes 

tables,  respective?;  ® ® 3/4\  and  7 3/8"  billet 
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7.  A bad  bay  priority  scheme,  I.e.,  the  bays  which  are  used  the  most  are 
located  at  the  far  end  of  the  yard  from  the  work  areas,  was  then  established 
and  the  simulation  was  again  repeated  three  times.  The  variability  between 
runs  was  less  than  3%  and  the  results  Indicate  that  the  crane  would  be  Idle 
at  least  9%  of  the  time  with  a probability  of  .995+.  The  .99  probability 
lead  times  required  by  the  feed  tables  are  11.13  minutes,  12.97  minutes, 

and  14.88  minutes  for  the  5 1/4",  6 3/4",  and  7 3/8"  billet  tables  respectively. 

8.  Chamberlain  objected  to  the  distributions  used  for  the  time  between 
charges,  because  they  did  not  allow  for  peaks  of  activity.  Therefore,  one 
run  was  made  with  these  distributions  changed  to  consist  of  a peak  at 'the 
mean  time  between  charges  when  the  lines  were  operating  at  340  mults-per- 
hour  and  a tail  to  the  right  to  account  for  downtime.  The  results  indicated 
that  the  idle  time  of  the  crane  dropped  from  33%  to  31%.  A single  run  was. 
judged  to  be  sufficient  since  the  previous  multiple  runs  had  shown  less  tfian 
3%  variation  in  the  results. 

9.  Additional  sensitivity  analysis  was  per*formed  by  varying  the  maximum 
speed  of  the  crane.  These  runs  were  made  with  both  good  and  bad  bay  pri- 
ority schemes,  normally  distributed  time  between  charges,  and  constant 
acceleration  at  1 foot  per  second  squared.  The  maximum  velocity  was  varied 
from  80  feet  per  minute  to  400  feet  per  minute,  and  the  results  are  plotted 
in  Figure  1.  If  a good  bay  priority  scheme  is  used,  then  top  speeds  in 
excess  of  200  feet-per-minute  only  Increase  the  reserve  capability  from 
33%  to  41%.  Figure  2 is  a plot  of  the  average  speed  of  the  crane  as  a per- 
centage of  its  maximum  speed.  With  a good  bay  priority  scheme,  the  maxi- 
mum usage  of  the  crane's  capabilities  occurs  at  approximately  160  ft/minute 
and  with  a bad  bay  priority  scheme,  the  maximum  occurs  around  240  ft/minute. 
Outside  of  this  range,  the  percent  utilization  of  the  crane's  speed  capabil- 
ity drops  off. 

10.  The  third  case  reported  in  the  original  report,  i.e.,  shock-loading 
the  rail  car  queue  with  a month's  supply  of  billets  and  determining  the 
time  required  to  unload  it,  was  not  performed  since  this  is  clearly  a situa- 
tion in  which  both  cranes  would  be  used  and  the  computer  model  does  not  have 
provisions  for  simultaneous  operation  of  the  two  cranes. 

DAVID  HOEHN 

Operations  Research  Analyst 

Logistics  Systems  Analysis  Division 
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1.  Reference; 

a.  MFR,  DRSAR-MA  to  DRSAR-SA,  9 Feb  76,  subject:  VADS  PIP  Provisioning. 

b.  FORECONs  between  DRSAR-MM,  Mr.  Crouch,  DRSAR-MA,  Mr.  Stehn,  and  DRSAR-SA, 

CPT  Krueger,  Mar  76  - Dec  76,  subject  as  above. 

c.  MFR,  DRSAR-MM,  10  Mar  76,  subject:  Initial  Provisioning. 

2.  The  Systems  Analysis  Directorate  was  tasked  (ref  la)  to  analyze  the  provisioning 
systeiTi  as  applied  at  IIQ,  ARRCOM.  The  attached  MFR  (Incl  1)  contains  the  analysis 

of  the  provisioning  system.  An  estimate  of  the  provisioning  time  was  developed  in 
March  1976  (ref  lb)  and  used  to  validate  the  developed  simulation.  A sutamary  of  the 
various  time  frames  simulated  is  shown  below: 


SIMULATIOri 

DAYS 

DIFFEREHCE 

a* 

Current  Process 

265 

- 

b. 

Increase  Provisioning  Personnel 

261 

-4 

c. 

Reduce  DISC  and  DSA/GSA  Time 

179 

-86 

d. 

ISA 

235 

-30 

e. 

Combinatorial  (bSd) 

230 

-30 

The  time  estimates  provided  by  provisioning  and  cataloging  personnel  (ref  lb  and  1c) 
were  used  as  input  for  the  simulation.  Even  though  the  above  changes  to  the  current 
process  show  reductions  in  processing  time,  they  should  not  be  implemented  until  cost 
savings  can  be  identified. 

3.  Point  of  contact  is  CPT  Krueger,  extension  6370. 


M.  RllIAN 

Director,  Systems  Analysis  Directorate 


1 Incl 
as 
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SUBJECT:  Analysis  of  the  Provisioning  System  ?s  Applied  at  HQ,  ARRCOM 


1.  Reference: 

a.  Minutes  of  C6,  ARRCOM  Weekly  Staff  Meet^ing,  /d  Jan  76, 

bi  DRSAR-MA  MFR  to  DRSAR-SA,  subject:  /vAD!;  PIP  Provisioning,  9 Feb  76 

' ' J / ■'  , ' ' i 

c.  Army  Regulation  700-18,  Provisioning  of 'US/Army  Equipment,  21  Sep/73 

I ■ ' L ' ' 

d.  Technical  Manual  38-715-1,  Provisioning  irechniques,  Oct  65. 

■ ■ I . If 

e.  Commodity  Command  Standard  System  Operating  Instructions,  Vol Vl, 

No.  18-700-13,  Provisioning  System,  Oct  75 J 

f.  Standard  Operating  Procedure  No.  700-MA-26,  Commodity  Corrmai 
Standard  System  Provisioning  System,  9 0ec'i74. 

g.  DRSAR-MM  MFR,  subject:  Initial  Provisioning,  10  Mar  76. 

, h.  Numerous  FONECONs  between  Mr,  Crojch,  DRSAR-MM,  Mr.  Stehn,  DRSAR-I^A, 
and  CPT  Krueger,  DRSAR-SAA,  subject:  Ref /la,  above,  period  of  time  - 
Mar  76  through  Dec  76.  ■ , / 

'/ . i 

i.  SAO  Note  2,  "Secondary  Items  Administrative  Lead  Time  Simulation 
Study,  Mr.  R.i  Banash,  etal , June  74. 

, 1 / i 

ij.  Military  Standard  No.  1388-1,  Logistical  Support  Analysis,  15  Oct  73, 

2.  Introdu 


ition. 


This  Directorate  was  tasked  initially  {ref  la)  to 


theiVULCAN  Air  Defense  System  (VADS)  Product  Improvement  Program 
(PIP)  provisioning  effort.  This  tasking  was  subsequently  broadened  to 

ArnPAu  overall  provisioning  system  as  applied  at  HQ, 

ARRCiOM,  (ref  lb).  j 

3.  Background.  Provisioning  is  the  process  for  determining  and  acquiring 
the  range  and  quantity  of  support  items  (repair  parts,  special  tools, 
technical  manuals,  etc)  necessary  to  operate  and  maintain  a weapon  system 
for  an  initial  period  of  service.  Provisioning  planning  begins  early  in 
the  life  cycle  of  a new  system  or  early  in  the  design  phase  for  product 
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improvements.  Hov'/ever,  adequate  guide! ines;  for  starting  this  planning 

and  establishing  the  sensitivity  of  the  provisioning  process  to  procedural 
factors  in  the  system  are  unknown.  The  purpose  of  this  study  was  to /analyze 
and  simulate  the  provisioning  system  as  prescribed  in  references  Ic^-lg 
and  perform  a sensitivity  analysis  on  identified  ^ystem  factors.  ' 

4.  Methodology . The  procedures  for  accomplishing  provisioning  actions  ' 
were  described  in  a procedural  flow  format  utilizing  major  contributions  ' 
from  the  Maintenance  and  Materiel  Management  Directorates  (ref  Ig  and  Ih). 
The  network  developed  related  activities  and  decision  points  within  the  ' 
provisioning  process.  This  network  was  then  simulated  using  the  General 
Purpose  Simulation  System  (GPSS)  developed  by  the  Science  Research  Asso- 
ciates, Inc.,, a subsidiary  of  IBM.  The  purpose  of  this  step  was  to 
obtain  an  automated  representation  of  the  provisioning  process  in  order 
to  quantitatively  assess  the  current  process.  After  verification  of  the 
process  simulation,  proposed  changes  and  identified  problem  areas  were 
analyzed.  A,ll  results  have  been  rounded  to  whole  days,  and  the  days 
refer  to  calendar  days. 

. ,1  ' ' _ 

Discussion  of  Data.  No  historical  data  existed  for  times  to  complete 
a provisioning  process;  however,  a generalized  time  estimate  was  developed 
in  March  1976  by  the  provisioning  and  cataloging  personnel  of  the  Mainte- 
nance and  Materiel  Management  Directorates  (ref  Ig).  This  estimate  was 
later  refined  (ref  Ih).  This  refined  estimate  was  used  to  verify  the 
simulation  results.  It  should  be  noted  that  this  study  covered  the 
provisioning  process  from  the  point  in  time  that  the  Maintenance  Direc- 
tordtG  TGceives  thG  initial  provisioning  input  from  the  developer  until 
that  point  in  time  that  all  National  Stock  Numbers  (NSN)  are  assigned. 

A separate  study  (ref  li)  simulated  the  administrative  lead  time  (ALT) 
for  the  procurement  of  secondary  items. 


^uj^>"ent  Process  Analysis.  The  current  process  was  simulated  using 
the  GPSS  computer  program.  GPSS  utilizes  three  basic  entities:  Facilities 
Storages,  and  Queues.  A facility  is  an  entity  that  can  handle  only  one 
transaction  at  a time,  for  example,  the  Configuration  Control  Board 
reviews  only  one  Engineering  Change  Proposal  at  a time.  A storage  entity 
IS  one  that  can  handle  up  to  and  including  a specified  number  of  trans- 
actions  at  one  time,  for  example,  the  provisioning  branch  can  handle  as 
many  transactions  as  it  has  personnel  available.  The  final  entity,  Queues, 

IS  an  entity  in  which  a transaction  waiting  to  enter  a facility  or  storage  - 
resides  until  space  is  available  for  it  to  be  processed. 


+ h three  entities  were  combined  in  a logical  manner  that  described 

the  actual  flow  process  for  provisioning  transactions.  The  system  simu- 
lation  was  allowed  to  reach  a steady  state  to  reflect  the  average  number 

rpcNif^  ^9  ^°'”P^®^®.3.^'"3nsaction.  Table  1 compares  the  simulation 
results  with  the  original  March  76  estimate. 
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■ TABLE  1 


Current  Process 


Simulation 


Estimate  (March  1976) 

I 

Difference 


65  (Calendar  Days) 


'I  ■ 


■ / / 


I ' 'I 

The  26  day  difference  between  the  simulation  resjult  and  the  original  '/ 
estimate  was  discussed  with  the  concerned  representatives  (ref  Ih),  Iv 
was  concluded  that  the  March  76  estimate  did  not' fully  allow  for  the  time 
needed  to  correct  errors  noted  when  the  Provisioning  Master  Data  Record 
was  received  from  the  ALPHA  system  and  reviev/ed*  ' Also,  the  possible're- 
submission  to  the  Defense  Logistics  Service  Center  (DLSC)  and  Defense  Supply 
Agency/General  Services  Administration  DSA/GSA  fpr  NSN's  was  not  fully 
taken  into  account.  It  was  felt,  therefore,  that|  the  simulation  run  was  ^ 
more  viable  time  estimate  and  would  serve  ^s  the  standard  against  which 
the  sensitivity  analyses  would  be  compared. 


■I 


/ 

7.  Sensitivity  Analyses.  The  sensitivity  of  the  following  system  factors 
to  change  was  investigated.  ' ' 


/ 


a.  Six  uhfilled  positions  in  the  Provisioning  Branch  of  the  Mainte- 
nance Directorate  had  been  identified.  The  opinion  was  that  this  shortage 
created  a backlog  situation.  | 

( i i 

, lb.  A reduction  of  the  minimum  required  process  time  at  DLSC  (DSA/GSA) 
from  60  to  |0  days  was  investigated. 

ic.  Excldsive  use 'of  LSA,  Logistical  Support  Analysis  was  investigated. 


d.  A combination  of  using  LSA  and  increasing  the  number  of  provisioning 
to  personnel  by  sik  was  investigated.  i 

No  increase  in  the  number  of  personnel  in  the  Cataloging  Division  was  con- 
sidered since  no  unfilled  positions  were  indicated  and  no  queue  time  develop* 
ed  in  the  simulation  that  would  indicate  the  possible  need  to  increase 
personnel . i 

!■  ' ' ~ 

8.  Sensitivity  Results. 


a.  Increase  of  the  Number  of  Personnel  in  the  Provisioning  Branch 

ioning  Branch  reported  six  unfilled  slots  exist 


• eaoc  vj  1 uiic  I 

The  Chief  of  the  Provisi( 


sted 
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within  his  organization.  This  fact,  coupled  with  existence  of  Queue  time 
for  his  organization, indicated  that  an  increase  in  personnel  may  have 
a potential,  to  decrease  the  provisioning  time  frame.  The  number  of 
personnel  for  the  Provisioning  Branch  was  increased  by  six  for  this 
simulation  run  and  no  other  changes  were  made  to  the  simulation  process. 

The  results  of  the  run  are  shown  in  Table  2. 

I ' ■ ' 

; ' TABLE  2 

!*  • . 

Increased  Number  of  Provisioning  Personnel 

Current  Process  . 265  (Calendar  Days) 

Increased  Number  261 

Difference  4 

f 

I 

While  the  addition  of  six  personnel  to  the  Provisioning  Branch  did 
decrease  the  time  frame  by  four  days,  it  is  doubtful  that  the  additional 
personnel  would  be  warranted.  To  further  substantiate  this,  the  increased 
number  of  personnel  was  varied  from  six  down  to  two  and  the  reduction  in 
provisioning  time  was  only  changed  from  four  days  to  three.  The  addition 
of  only  one  person  resulted  in  a reduction  of  approximately  one  day. 

b.  Reduction  in  DLSC  and  DSA,/GSA  Process  Time.  It  was  pointed  out 
that  DLSC  and  DSA/GSA  are  allowed  60  days  to  process  requests  for  stock 
numbers.  This  is  virtually  a fixed  time  delay  and  if  some  stock  number 
data  were  found  to  be  in  error,  an  additional  60  days  would  be  allowed 
DLSC  or  DSA/GSA  upon  resubmission  of  the  data.  While  a reduction  in  this 
time  frame  is  not  within  HQ,  ARRCOM' s control , the  results  shown  in  Table  3 
may  assist  in  causing  revision  of  this  60  day  allowance.  The  reduction 
used  in  this  simulation  was  a cut  of  one-half,  60  to  30  days,  and  all  other 
\ portions  of  the  simulation  remained  unchanged. 

TABLE  3 

Reduced  Process  Time  DLSC  and  DSA/GSA 
Current  Process  265  CCalendar  Days) . 

Reduced  DLSC  and  DSA/GSA  Time  179 
Difference  '86 

The  reduced  times  made  a significant  contribution  to  reducing  the  time 
required  to  properly  provision  a weapon  system.  If  the  60  day  limit  is 
absolute,  time  still  could  be  saved  in  provisioning  if  error  resubmission 
could  be  placed  in  a category  that  waived  the  mandatory  60  day  process  time. 
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c.  Logistical  Support  Analysis  (ISA).  A Logistical  Support  Analysis 
(LSA)  was  to  be  required  on  all  newly  developed  weapon  systems  and  major 
modifications.  This  requirement,  however,  is  not  being  fully  implemented 
because  of  some  difficulties  in  the  bridging  program  that  converts  the 
LSA  input  to  the  CCSS  format,  and  because  of  the  need  to  thoroughly  familiar- 
ize personnel  with  the  use  of  LSA,  t.e.,  a developer  (military  and  civilian), 
personnel  and  provisioning  and  cataloging  personnel.  The  procedures  to 
use  LSA  In  the  provisioning  process  were  simulated  and  the  results  are  , 
shown  in  Table  4. 


TABLE  4 


Use  of  Logistical 

Support  Analysis 

Current  Process 

' 265  (Calendar  Days) 

LSA  Process 

235 

Difference 

30 

The  use  of  LSA  is  completely  within  ARRCOM’s  control  to  implement  and 
training  in  the  use  of  LSA  is  available.  If  cost  savings  can  be  identi- 
fied with  this  30  day  reduction,  then  this  alternative  may  warrant  imple 
mentation. 


d.  Combinatorial . The  use  of  LSA  and  an  increase  of  six  in  the  current 
number  of  provisioning  personnel  was  simulated  next.  This  combination  v/as 
chosen  because  of  the  ability  of  HQ,  ARRCOM  to  readily  implement  these 


changes. ; 

1 

. 

i 

. 1 ' 

Current  Process 

Combinatorial 

Difference 

TABLE  5 

Combinatorial 

!!■ 

■ i 

■ r 

1 

j.|  .. 

265  '(Calendar  Days) 
235!  'i  ' 

30  • 

The  results  of  this  simulation  run  are  identical  with  just  using  LSA. 
Further  investigation  revealed  that  the  Provisioning  Branch  realized  no 
queue  time  under  just  LSA  procedures  and  the  processes  needed  to  conduct 
a provisioning  effort  still  require  the  same  amount  of  time.  Thus, 
increasing  the  available  number  of  personnel  when  no  backlog  exists  would 
riot  create  a reduction  in  time. 

9.  Summary.  The  above  results  indicate  that  the  least  effective  means 
of  reducing  provisioning  time  frames  is  to  fill  all  or  a portion  of  the 
unfilled  positions  in  the  Provisioning  Branch.  The  most  effective  means 
of  reducing  the  provisioning  time  is  to  have  DLSC  and  DSA/GSA's  mandatory 
\ time  frame  reduced;  however,  this  is  not  an  immediate  solution,  since 

••  I • . » 
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it  is  an  action  required  at  DOD  level.  The  combinatorial  alternative  is 
not  a reasonable  solution  since  the  addition  of  personnel  resulted  in  no 
change  over  using  ISA  alone.  The  most  viable  solution  is  to  use  LSA.  To 
do  this  may  require  special  effort  to  insure  the  bridging  program  problems 
are  solved  and  that  the  largest  possible  number  of  personnel  are  trained 
to  use  LSA. 
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M551  Sheridan  Wear-Out 
DRSAR-SA 


DRSAR-HA 

DR5AR-M.M 


6 Jl^  ^977 

CPT  Krueger/jl s/6370 


1.  Reference: 

a.  DF  w/incl,  DRSAR-MA  to  DRSAR-SA,  1 Nov  76,  subject:  Trip  Report  - USAREUR 
Visit  11-25  Oct  76,  DCG,  Dir  of  Maintenance  and  Ch,  DRSAR-MMP,  Part  II. 

b.  Demand  Return  Disposal  Files,  15  Dec  74-15  Dec  76. 

c.  Order  of  Merit  List,  11551  Sheridan,  6 Oct  76. 

2.  The  Systems  Analysis  Directorate  was  tasked  (ref  la)  to  conduct  a failure  rate 
analysis  on  European  based  Sheridans  and  additionally,  to  determine  if  the  component 
buy  policy  for  ARRCOM  managed  M551  Sheridan  components  and  repair  parts  need  to  be 
revised  to  accommodate  a rapidly  increasinvg  wear-out  rate.  As  discussed  in  the 
attached  IlFR  (Incl  1 ) , no  rapidly  increasing  wear-out  rate  was  indicated.  As  a result, 
this  study  indicates  that  no  alteration  of  the  component  buy  policy  is  needed. 

3.  In  regard  to  the  processes  for  collecting  data  on  parts  demand  history,  tvra 

important  changes  need  to  be  made.  First,  the  Demand  Return  Disposal  Files  (ref  lb) 
only  contain  the  past  two  years  demand  history.  The  number  of  years  of  demand 
history  retained  at  HQ,  ARRCOM  on  magnetic  tapes  needs  to  be  extended  to  a minimuni 
of  at  least  four  years.  Second,  the  weapons  codes  (per  AR  725-50)  used  in  the  DRD\ 
files  are  not  used  in  Europe  and  are  not  a required  entry  of  CONUS.  The  lack  of 
these  codes  hinders  the  uae  of  the  DRD  files  as  an  accurate  data  base.  The  feasi- 
bility of  making  the  use  of  those  weapons  codes  prescribed  in  AR  725-50  a mandatory  \ 

Army  wide  entry  should  be  investigated. 

4.  Point  of  contact  is  CPT  Krueger,  extension  6370. 


as 


I Incl 


M.  RHIAN 

Director,  Systems  Analysis  Directorate 
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1.  Reference: 

a.  DRSAR-MA  DF  w/incl  to  DRSAR-SA,  1 Nov  76,  subject:  Trip  Report  - 
USAREUR  Visit  11-25  Oct  76,  DCG,  Dir  of  Maintenance  and  Ch,  DRSAR-Mf>iP, 
Part  II, 

b.  Demand  Return  Disposal  Files,  15  Dec  74-15  Dec  76.  ■ 

I 

c.  Order  of  Merit  List,  M551  Sheridan,  6 Oct  76.. 

2.  Introduction.  The  Systems  Analysis  Directorate  was  tasked  to  conduct 
a failure  rate  analysis  on  European-based  Sheridans  and  additionally,  to 
determine  if  the  component  buy  policy  for  ARRCOM  managed  Sheridan  related 
items  need  be  revised  to  reflect  an  accelerated  failure  rate  for  Sheridan 
turret  components  and  repair  parts. 

3.  Background.  During  a visit  to  USAREUR  (ref  la),  it  was  reported  to 
the  DCG,  ARMCOM  that  the  age  of  the  Sheridans  in  USAREUR  was  a factor 
causing  a rapidly  increasing  component  wear-out  rate.  This  study  was 
approached  from  the  viewpoint  that  if  the  European-based  Sheridans,  all 
having  varied  ages,  were,  in  fact,  v/earing  out'at  a rapid  rate,  the  com- 
ponent buy  policy  should  be  adjusted  accordingly. 
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1 Vo// 


DRSAR-SAA 

SUBJECT;  M551  Sheridan  Wear-Out- 


5.  Methodology.  The  rationale  for  using  the  DRD  files  vms  that  if  the 
variably-aged  Sheridans  were  rapidly  wearing  out,  this  would  be  reflected 
in  a corresponding  rapid  increase  in  the  number  of  demands  and/or  quanti- 
ties demanded.  A COBOL  program  was  written  that  would  first  screen  the 
DRD  files  for  the  stock  numbers  obtained  from  the  OML  and  then  further 
screen  the  remaining  stock  numbers  for  recurring  European-generated  demands. 
This  last  screening  process  reduced  the  stock  numbers  to  be  analyzed  to 
as  close  to  the  using  unit  level  as  possible. 


The  resultant  European-identified  Sheridan  , stock  numbers  were  then 
subjected  to  four  analyses  to  determine  if  a rapid  increase  had  occurred 
in  the  number  of  demands  and/or  quantity  demanded. 


a.  The  stock  numbers  were  first  ordered  by  .total  cost  per  NSN  from 
high  to  low.  The  top  90%  of  the  cumulative  total  cost  were  then  subjected 
to  a trend  analysis.  This  trend  analysis  consisted-  of  comparing  first 
and  second  year  demands. 


b.  The  total  demands  per  month  and  total  quantity  demanded  per  month 

were  each  graphed  (Incl  1 A 2)  and  again  a trend-,  analysis  performed  on  the/ 
data.  , i I ■ 

I ’ • 1 1 . 

c.  The  average  quantity  per  demand  was  graphed  (Incl  3)  by  month  an/ 
again  a trend  analysis  was • performed  on  the  data: 

j 

d.  The  average  cumulative  quantity  demanded  'per  month  was  calculated, 

graphed  (Incl  4)  and  a trend  analysis  performed.  - 

6.  Results  of  Analyses.  ' '■ 


a.  After  the  European  stock  numbers  v/ere  identified  and  subsequently 
ranked  by  cost  from  high  to  low,  the  stock  numbers  comprising  the  top  90% 
of  the  total  cost  were  separated  for  analysis.  The  number  of  NSN's  in  this 
category  was  111  or  11.6%  of  all  the  stock  numbers  considered.  In  addition 
these  stock  numbers  accounted  for  35.4%  of  the  total  demands  for  all  stock 
numbers  considered.  The  number  of  demands  for  the  first  and  second  years 
for  each  stock  number  was  determined  and  a linear  regression  performed 
in  order  to  determine  the  trend  from  the  first  year  to  the  second  year. 

The  linear  regression  yielded  a positive  slope  which  would  indicate  an 
increase  from  one  year  to  the- other  and  the  correlation  factor  (a  measure 
of  how  close^  the  data  points  conform  to  a straight  line)  was  moderately 
high,  v/hich  'indicates  a close  fit.  However,  it  was  felt  that  further 
investigation. was  warranted.  This  need  became  more  obvious  when  a closer 
examination  of  the  data  was  made.  First,  the  demands  were  grouped  by  years 
as  opposed  to  months,  in  addition,  there  were  many  anamolies  in  the 
data,  i.e.,  demands  greatly  increased  or  decreased  from  the  first  to  the 
second  year.  Second,  upon  examination  of  specific  stock  numbers  that 
exhibited  a large  increase  in  demands  for  the  second  year,  it  was  found 
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that  this  increase  was  due  a great  extent  to  Maintenance  Letters  that  had 
been  sent  to  the  field  directing  various  maintenance  or  inspection  actions. 
These  letters  invariably  lead  to  a temporary  increase  in  demands.  In  addi- 
tion, almost  all  the  increased  second  year  demand  stock  numbers  v/ere' 
included  in  pending  Engineering  Change  Proposals  (ECPs)  for  the  Sheridan. 

In  all,  no  conclusions  could  be  made  as  to  any  indication  of  rapid  wear- 
out. 

' I 

b.  ^ The  total  demands  and  total  quantity  demanded  were  graphed  by 
month  in  an  attempt  to  determine  where  and  when  the  above  noted  fluctua-  ■ 
tions  occurred  and  also  to  gain  a better  insight  into  the  nature  of  the 
demands  i.e.,  ms  there  a rapid  increase  in  demands.  Because  this  data 
was  a time  series  analysis,  the  series  was  decomposed  to  eliminate  many 
of  the  anamolies  described  above,  such  as  the  seasonal  and  cyclical 
fluctuations  in  order  to  arrive  at  the  secular  trend  (the  secular  trend 
indicates  the  long-run  growth  or  decline  of  the  series).  Upon  analyzing 
the  secular  trend,  no  indication  of  rapid  increase  v/as  detected. 

c.  The  average  quantity  per  demand  v/as  plotted  by  month  and  again 

a secular  trend  analysis  was  made  to  determine  if  a rapid  rise  was  indi- 
cated. This  attempt  again  failed  to  indicate  anything  other  than  a 
slight  increase. 

d.  The  average  cumulative  quantity  demanded  per  month  was  calculated 
and  graphed.  This  was  another  means  of  trying  to  discern  a rapid  increase 
in  demand  for  parts.  Upon  reviewing  the  results,  again  there  was  no 
indication  of  rapid  rise  in  demands  or  quantity  demanded. 

7.  Surnmary. 


a.  All  attempts  to  demonstrate  that  the  European-based  Sheridans  v;ere 
wearing  out  rapidly  failed  to  indicate  any  such  pattern.  A gradual  increase 
was  detected;  however,  this  increase  is  basically  inconclusive  because 

of  the  availability  of  only  two  years  of  demand  data  to  analyze.  The 
largest  benefit  to  be  realized  from  this  effort  is  that  the  DRD  files  are 
a source  of  excellent  data  provided  the  following  two  conditions  are  met. 
First,  more  than  just  two  continuous  years  of  demand  history  is  needed 
and  second,  that  all  weapons  codes,  as  assigned  in  AR  725-50,  be  a required 
entry  on  all  requisitions  originating  in  CONUS  and  from  OCONUS. 

/ i 

b.  From  the  analysis  above  there  appears  to  be (based  on  available 
data  ) no  need  to  alter  the  current  component  buy  policy.  However,  with 
at  least  four  years  of  demand  data,  the  above  analyses  may  be  able  to 
adequately  predict  when  changes  may  be  needed  before  problems  arise. 
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QUANTITY  DEMANDED/NUMBER  OF  REQUISITIONS 


Ciunulatlve  Average  Monthly 
Quantity  Requested 


I 


120 

115 

no 

105 

100 

95 

90 

85 

80 

75 

70 

65 

60 

55 

50 

^5 

40 

35 

30 

25 

20 

15 

■ 5 . 


4 


4 


i 


1D8 


DRSAR-SAM 


* 3 JUN  m 


MEMORANDUM  FOR  RECORD 

SUBJECT:  Statistical  Methods  Pertinent  to  a Potential  Ignition  Problem  in 
the  M188E1  Propelling  Charge 


1.  References: 

a.  MFR,  STEAP-MT-G,  11  Apr  77,  subject:  M188E1  Charge— Max  Negative 
aP  Results. 

b.  Letter,  DRDAR-LCU-E-P,  24  May  77,  subject:  Test  Program  Request 
ADEP  2061  Charge  Propelling  8-Inch  M188E1. 

2.  Background 

The  author  was  asked  by  the  PM-110E2  to  review  a potential  safety  prob- 
lem in  the  MllOAl  SP  howitzer  when  using  the  M188E1  propelling  charge. 
During  development  testing  of  this  charge,  pressure  is  measured  simultane- 
ously at  two  locations:  near  the  breechface  and  at  the  forward  end  of  the 
chamber  (tov/ard  the  muzzle).  After  proper  ignition  and  throughout  the 
interior  ballistic  cycle,  the  pressure  is  generally  higher  at  the  breech 
than  at  the  base  of  projectile.  However,  during  ignition  reverse  pressure 
gradients,  i.e.,  a larger  pressure  forward,  can  occur.  The  magnitude  of 
the  negative  pressure  differential,  measured  in  the  manner  indicated  above, 
has  been  found  to  correlate  positively  with  the  peak  chamber  pressure  sub- 
sequently experienced  during  the  interior  ballistic  cycle.  Thus,  a large 
absolute  pressure  differential,  Ap,  is  accompanied  by  a large  value  of 
the  peak  chamber  pressure. 

3.  Statement  of  the  Problem 


Because  pmax  limited  to  a value  consistent  with  projectile  and 

cannon  allowable  stresses,  it  follows  that  the  associated  value  of  Ap  must 
be  limited  by  some  safe  value.  The  Ap  limit  is  determined  in  part  by  the 
relationship  between  p^ax  Ap.  Presently  this  relationship  is  poorly  de- 
fined and  may,  in  fact,  oepend  upon  other  variables  such  as  propellant 
temperature.  Even  if  the  safe  limit  of  Ap  were  well  defined,  the  risk  of 
exceeding  this  limit  depends  upon  the  probability  distribution  function  of 
the  random  variable  Ap,  which  itself  may  depend  upon  propellant  temperature, 
among  other  variables.  Data  analyzed  in  Ref  a.  indicate  that  the  probability 
distribution  of  Ap  is  affected  by  large  changes  in  propellant  temperature. 
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During  safety  tests  of  the  MHO  system,  the  M188E1  propelling  charge  had 
been  subjected  to  a sequence  of  rough  handling  operations  and  subsequently 
temperature  conditioned  to  -50  degrees  F,  Several  of  the  rounds  fired  with 
this  treatment  experienced  Ap  values  in  excess  of  2 ksi  and  displayed 
values  of  p^^^  which  generally  increased  with  Ap,  Furthermore,  tests  of 
M188E1  charges  with  intentionally  mal constructed  igniters  when  subjected  to 
-50  degree  F conditioning  have  displayed  extremely  large  (>5  ksi)  values  of 
Ap  and  anomolously  large  values  of  p 

5,  In  the  light  of  these  results  STEAP-MT  has  begun  a series  of  tests 
(Ref.  b)  whose  general  purpose  is  to  better  quantify  the  factors  to  be  con- 
sidered in  estimating  the  risk  of  a catastrophic  malfunction  of  the  M188E1 
propelling  charge.  Since  the  anomolously  large  chamber  pressures  accompany 
improper  ignition  at  low  temperature,  it  is  important  to  define  the  effect  of 
temperature  on  the  probability  distribution  of  Ap  in  unmodified  charges 

lowing  anticipated  operational  rough  handling, 

6,  A Particular  Issue 


In  this  connection  an  immediate  question  is  whether  the  probability 
distribution  of  Ap  has  a temperature  dependence  for  temperatures  below  zero 
degrees  F,  To  answer  the  question  of  temperature  dependence  as  efficiently 
as  possible,  one  requires  powerful  statistical  tests  and  should,  of  course., 
make  use  of  all  applicable  existing  data.  With  these  things  in  mind,  I 
have  prepared  some  statistical  methods  which  may  be  helpful  in: 

a.  selecting  a statistical  sample  to  provide  an  adequate  degree  of 
discrimination  between  propelling  charge  treatments, 

b,  estimating  the  value  of  Ap  which  would  be  exceeded  at  a given  risk 
(a  percentile  of  the  distribution  of  Ap)  and  an  associated  confidence  inter- 
val for  this  estimate. 

7,  The  derivation  of  some  pertinent  statistical  tests  and  the  presentation 
of  their  operating  characteristics  for  several  sample  sizes  are  given  in 
Attachment  1 (Incl  1),  The  treatment  is  not  intended  to  be  exhaustive  but 
rather  to  define  the  power  of  some  parametric  and  non-pa rametric  tests  for 
this  particular  application.  Computer  programs  are  presented  in  Attachment  2 
(Incl  2). 
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SOME  STATISTICAL  METHODS  PERTAINING  TO  THE 
IGNITION  PROBLEM  IN  THE  M188E1  PROPELLING  CHARGE 


Background 

At  lov/  temperatures  following  sequential  rough  handling  tests  of  the 
XM188  propelling  charge,  some  unusually  large  maximum  absolute  values  of 
negative  pressure  differential  have  been  experienced  within  the  combustion 
chamber  of  the  eight-inch  M201  cannon.  Typically,  a large  value  of  the 
maximum  absolute  pressure  differential,  Ap,  accompanies  a large  peak 
chamber  pressure  during  the  interior  ballistic  period.  Several  tests 
have  been  proposed  to  investigate  this  problem,  which  is  regarded  as 
potentially  serious  because  of  possibly  unsafe  peak  chamber  pressures. 

In  one  series  of  tests  it  is  proposed  to  fire  sets  of  these  charges 
conditioned  at  several  temperatures,  in  order  to  determine  the  effect  of 
temperature  on  the  probability  distribution  of  Ap  produced  during  igni- 
tion of  the  charge. 

Parametric  Methods 

Available  data  on  the  XM188  charge  indicate  that  at  both  low  (-50°F) 
temperatures  and  high  temperatures  (145°F)  a two-parameter  Weibull  distribution 
is  a reasonable  statistical  model  for  Apt  These  data  suggest  that  the 
shape  parameter,  b,  is  nearly  unity--the  high- temperature  distribution 
having  a value  of  B only  slightly  in  excess  of  1,  and  the  low-tempera- 
ture distribution  having  a B slightly  less  than  1.  However,  due  to  the 
limited  sample  these  distributions  are  not  statistically  distinguishable 
from  the  (negative)  exponential  distribution  which  corresponds  to  a 
Weibull  with  B = 1.  Therefore,  in  the  following  analysis  the  exponential 
model  is  assumed.  Because  of  its  relationship  to  risk,  the  first  topic 
addressed  is  percentile  estimation. 

Variance  of  Percentile  Estimates  for  Exponential  Random  Variables 

The  lOOp  th  percentile,  Xp,  of  an  exponential  distribution  with  para- 
meter e is  given  by 

*MFR,  STEAP-MT-G,  11  Apr  77,  subject:  M188E1  Charge --Max  Neg  AP  Results. 
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or 


p = 1 - exp  - Xp/0 


Xp  = -0  ln(l  - p) 


(1.1) 

(1.2) 


Thus  with  p chosen,  x is  proportional  to  the  parameter  0. 

P A . 

The  maximum  likelihood  estimate  of  0,  0,  obtained  from  a sample  of 
n:  {X^. , i = 1,  n}  is,  simply, 

0 = i:",!  X./n  . (1.3) 

To  obtain  an  estimate  of  the  variance  of  x , one  notes  that,  from 

(1.2), 


Var(Xp)  = Var(0)  ln^(l  - p)  . 


(1.4) 


Now,  a well  known  result  for  the  exponential  distribution  is  that  the 
parameter  2n0/0  has  a chi-squared  (X-  ) distribution  with  2n  degrees  of 
freedom.  Further, 


Var(X-^^)  = 4n  , 


(1.5) 


so  that 


Var(0)  = 0vn  . 

Thus,  from  (1.4)  and  (1.6)  , 

Var(Xp)  = 0^  In^d  - p)/n  (1.7) 

or,  with  (1.2), 

Var(Xp)  = Xp/n  . (1.8) 

^ 2 

Note  that  the  distribution  of  2nx  /x  is  also  X-  with  2n  degrees  of 

P P 

freedom  since  Xp  and  0 are  proportional. 


Detecting  a Difference  Between  Two  Samples 

Suppose  that  two  samples  of  n each  are  used  to  estimate  the  x th 

A A P 

percentile  (and  parameter  0),  producing  0 and  0^.  Making  use  of  the 
Central  Limit  Theorem,  for  n greater  than  about  10  the  following  statis- 
tic is  approximately  standard  normal: 


6 - Cq  - (e  - Op) 
[Var(6)  + Var(ep)]''^^ 


(1.9) 
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Alternatively,  with 


f = 6/e^  , 

(1.10) 

(e  - e^)*^  (f  - 1)*^ 

(1.11) 

[e^  + 00]^"^^  (f^  + 1)^"^^ 

One  hypothesizes  that  6 = 

With  the  gaussian  assumption. 

e.  , with  the  alternative  H, : e > e . 
0 1 0 

Pit  < 1.65}  = 0.95  . 

(1.12) 

Consequently,  at  a risk  of  only  5% 
can  accept  H-j  if 

of  the  declaring  false  if  true,  one 

(e  - eJ/fT 
> 1.65  . 

7^2  7^2 

*'6  +0^ 

0 

(1.13) 

Notationally,  let  e = P{accept 
Then, 

(1.14) 

(e  - eJ/F 

e = P{ 2 1.65  > 0} 

/ 0 + 

(1.16) 

Using  the  expected  value  of  the  denominator  in  (1.15),  approximately, 
e = PU  + - 1 .65  > 0} 

(r  + l)‘/2 


or 

e = $(/fT(f  - l)(f^  + l)"''/^  - 1.65),  (1.16a) 

with 

1 y2/2 

♦(z)  = — e dx  . (1.16b) 

Plots  of  e(f)  for  various  values  of  n,  calculated  from  (1.16),  are  dis- 
played in  Figure  1. 
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By  a Factor  (f)  Between  Two  Sets  of  Experiments  with  Sample  Size  in  Each 


These  operating  characteristics  or  OC  curves  indicate  the  value  of  the 
fractional  shift  f which  must  occur  to  detect  a shift  in  the  parameter 
e between  samples.  The  rule  used  to  detect  a shift  is  given  by  (1.13) 
with  e and  6^  calculated  using  equation  (1.3).  The  results  of  several 
experiments  can  be  compared  by  applying  the  above  test  to  all  possible 
pairs,  where  the  basic  parameter  6^  is  estimated  from  the  sample  whose 
population  parameter  e is  expected  to  be  minimal  on  physical  grounds. 
However,  if  the  number  of  sets  of  samples  at,  say,  different  temperatures 
is  large,  pairwise  comparison  is  not  the  most  powerful  statistical 
method  to  detect  a temperature  effect.  If  the  number  of  experiments  is 
greater  than  about  3,  regression  of  e on  temperature  appears  preferable. 

The  use  of  parametric  methods  to  discriminate  between  treatments 
requires  an  assumption  concerning  the  form  of  the  distribution  function. 
If  the  nature  of  the  distribution  function  is  seriously  in  doubt,  par- 
ticularly for  large  values  of  the  argument  Ap,  as  is  the  case  here,  it 
is  preferable  to  use  non-parametric  methods.  In  the  next  section  a 
specific  non-parametric  statistic  is  suggested  for  detecting  the  effect 
of  treatm.ent  when  three  treatments  are  applied  to  three  samples.  To 
facilitate  comparisons  between  the  operating  characteristics  of  the  above 
parametric  test  and  the  non-parametric  tests,  we  display  in  Figure  2 
the  relationship  between  f — as  defined  in  (1.10)--  and  a probability 
used  in  the  non-parametric  tests,  namely,  the  proportion  of  the  popula- 
tion of  Ap  lieing  above  2000  psi. 
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Non-Parametric  Methods 


To  detect  whether  the  experimental  treatment  affects  the  distribu- 
tion of  Ap  (the  variable  of  interest),  one  can  examine  some  non-para- 
metric  measure  such  as  the  number  of  occasions  in  a sample  of  n in  which 
Ap  exceeds  a specified  value.  This  reduces  the  analysis  of  each  experi- 
ment to  counting  the  number  of  such  "successes"  over  Bernoulli  trials. 

The  discrimination  between  treatments  is  based  upon  a comparison  of  the 
number  of  successes.  If  one  is  merely  interested  in  whether  any  sort  of 
change  occurs  in  the  probability  of  success  over  the  three  sets  of 
Bernoulli  trials,  one  possible  test  statistic  is  the  (discrete)  binomial 
range  or  extreme  spread*  in  the  number  of  successes.  This  statistic  has 
the  virtue  that  it  combines  the  results  of  all  experiments  into  one  index 
value.  In  the  following  derivations  the  distribution  function  for  the 
binomial  range  is  developed  and  used  to  formulate  a test  to  detect  depar- 
tures from  constancy  in  the  probability  of  success  over  the  three  sets  of 
Bernoulli  trials,  i.e.,  to  detect  ^3-  The  operating  character- 

istic of  this  test  is  calculated.  In  this  context  the  operating  charac- 
teristic is  the  probability  of  detecting  a departure  from  constancy  as  a 
function  of  the  magnitude  of  the  departure. 


The  Distribution  for  the  Discrete  Range  of  Successes  in  Three  Sets  of 
Bernoulli  Trials 

Given  the  event  E^(i):  during  the  i th  experiment,  there  are  r 
rounds  which  have  a Ap  in  excess  of,  say,  2 ksi,  given  n rounds  per  experi- 
ment are  fired  with  the  probability  equal  to  77.  that  for  any  round  of  the 
i th  set  Ap  will  exceed  2 ksi.  Then, 


P{E^(i)}  = A^(1)  = (J)Tr^(l  - 77.)" 


-r 


(2.1) 


0 r _<  n . 

The  expected  number  of  rounds  in  the  i th  experiement  with  "success," 
i.e.,  having  a Ap  in  excess  of  2 ksi,  is  n 77^.. 


weenX^UrgesranfsLnesrvIlJlfof  rsamjle!''  difference  bet- 
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The  range  or  maximum  spread  (s)  of  the  observed  successes  from  all 
of  the  three  experiments  has  a probability  density  (p.d.f.)  dependent 
upon  n:  P3(n),  with  s a member  of  the  discrete  set  S, 

S~{0jlj4..jn}.  (2.2) 

The  p.d.f.  for  the  range  is  derived  as  follows.  Let  r be  a dummy 
variable  for  the  number  of  successes  in  the  first  experiment.  Then, 

P{spread  is  exactly  s}  = p^  . 

The  expression  for  p^  is  developed  by  exhaustively  enumerating  the  events 
which  produce  a maximum  spread  s and  then  writing  the  probabilities  of 
these  events  and  taking  their  sum. 

Ps  ^ ^r+s^^^’  given  r+s  e S}- 

y{Eq(3)}  + P{E^_3(3)  or  E^^^(3),  given  r+s  e S}ZqP{Eq(2)} 
-P{E^_^(2)  or  E^^^(2),  given  r+s  e S}- 

where  the  limits  of  q are 

V-„  = max  (0,  r - s) 

"max  ' <"■  ^ ^ • 

In  evaluating  the  second  factor  within  the  sum  on  the  r.h.s.  of  (2.3) 
only  events  for  which  the  indices  are  in  the  set  S are  evaluated. 

Example:  n = 3 

Pq  = Aq(1)Aq(2)Aq(3)  + A.,(1)a.,(2)a.,(3)  + A2(1)A2(2)A2(3)  + 

X3(1)A3(2)A3(3)  (2.4) 

Pi  = Ao(1)A.j(2)Aq(3)  + Aq(1)A.|(3)(Aq(2)+X.|(2))  + 

A-,(1)(Aq(2)+A2(2))a.,(3)  + A.,(1)(Xq(3)+X2(3))(Aq(2)+A.,(2)+A2(2))  + 
A2(1)(a^(2)+A3(2))X2(3)  + A2(1)(A^(3)+A3(3))(a^(2)+A2(2)+X3(2))  + 
A3(1)X2(2)X3(3)  + X3(1)X2(3)(X2(2)+X3(2))  (2,5) 
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^i(1)A3(2)(a^(3)+A2(3))  + A^(1)A3(3)  + 

A3(1)A^(2)(A2(3)+A3(3))  + A3(1)a^(3)(a^(2)+A2(2)+A3(2))  (2.6) 

P3  = Aq(1)A3(2)  + Aq(1)A3(3)(Aq(2)+A^(2)+A2(2))  + 

A3(1)Ao(2)  + A3(1)Aq(3)(a^(2)+A2(2)+A3(2))  (2.7) 

A numerical  evaluation  of  (2.4)  thru  (2.7)  with  tt  = u = u = 0.1 
produces  the  following  distribution  of  binomial  range:  ^ 


s 

*^s 

"^1  = 1 Pi 

0 

0.4018 

0.4018 

1 

0.5315 

0.9332 

2 

0.0644 

0.9976 

3 

0.0024 

1 .0000 

Summary  statistics  are: 

mean  = sj^iip.  = 0.66735 

std.  dev.  = [lj^.|i^p^.  - mean  ^ 0.60418 

coefficient  of  variation  = std.  dev. /mean 

= 0.90535 
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Test  for  Constancy  of 

The  distribution  for  the  discrete  binomial  range  was  evaluated,  using 
(2.3),  for  several  sets  of  values  of  the  parameters  n,  tt-j  , Tr2,  These 
results  are  shown  in  Table  1.  One  notes  that  a progressive  departure 
from  constancy  of  the  tt's  in  the  manner  indicated  in  Table  1 is  accompanied 
by  a shift  in  the  distribution  of  the  binomial  range  to  the  right,  i.e., 
in  the  direction  of  larger  values.  Further,  even  though  the  standard  devi- 
ation also  increases  with  increasing  or  (tt^  being  fixed),  the  coeffi- 
cient of  variation  decreases.  Thus,  the  distribution  becomes  relatively 
less  disperse. 

These  characteristics  of  the  distribution  of  binomial  range,  s, 
suggest  a simple  test  of  constancy  of  the  tt's.  Specifically,  for  a given 
sample  size  (n),  select  a value  of  s for  which  the  tt's  would  be  declared 
identical.  Call  this  the  acceptance  number  a.  For  values  of  s greater 
than  a,  one  would  accept  the  alternative  hypothesis,  viz.,  the  tt's  are  not 
all  identical.  That  is,  if  s > a,  the  treatment  is  declared  to  affect 
the  value  of  tt,  the  probability  that  Ap  exceeds  2 ksi. 

In  selecting  the  acceptance  number  a for  a given  n,  one  must  decide 
what  risk  will  be  accepted  in  declaring  the  tt's  different  if  they  are  in 
fact  not.  For  example  for  n = 10  and  a = 2,  from  Table  1,  the  risk  is 
about  10%.  Similarly,  for  n = 20  and  a = 3 this  risk  is  about  11%,  and 
for  n = 30  and  a = 4 the  risk  is  approximately  9%. 

If  a sample  (n)  of  30  and  an  acceptance  number  (a)  of  4 were  chosen, 
the  probability  of  detecting  a shift  of  and  from  0.1  to  0.4  would 
exceed  98%  (from  Table  1).  It  is  of  interest  to  compare  the  power  of 
this  test  to  that  of  the  previous  parametric  test  on  the  difference 
6-00  . To  facilitate  the  comparison,  note  from  Figure  2 that  the  value 
of  f corresponding  to  an  ordinate  of  0.4  (=  n)  is  2.513.  Then,  using  this 
value  of  f and  assigning  the  same  risk  of  mistaking  a shift  of  tt  as  in  the 
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non-parametric  test,  viz.  0.0927,  the  parametric  test  would  yield  a 
96%  probability  of  detecting  this  shift.  In  this  case  the  non-para- 
metric test  is  actually  more  discriminating.  The  reason  for  the  better 
discrimination  of  the  non-parametric  test  is  that  it  takes  information 
from  all  three  experiments  rather  than  from  simply  a pair  as  does  the 
parametric  test. 

Operating  Characteristic  of  the  Non-Parametric  Test 

Using  the  results  in  Table  1 (with  equal  values  of  -n^  and  n^),  one 
can  develop  the  operating  characteristics  of  three  non-parametric  tests 
using  the  binomial  range  with  values  of  n = 10,  20,  30  and  corresponding 
values  of  a = 2,  3,  4.  The  probability  of  accepting  the  alternate 
hypothesis  (H-i)  that  the  tt's  are  different  is  shown  in  Figure  3 as  a 
function  of  ^2  and  tt^.  In  constructing  Figure  3,  it  is  assumed  that  the 
values  of  ^2  and  are  the  same.  In  this  case  the  probability  of  accept- 
ing has  only  a single  argument.  However,  this  assumption  is  somewhat 
restrictive.  In  general,  P{accept  depends  upon  two  arguments— it2  and 
which  may  be  different.  The  latter  relation  is  shown  in  Figure  4,  an 
isometric  graph.  It  is  noted  that  the  region  in  the  domain  of  ^2  and 
for  which  the  probability  of  accepting  is  less  than  0.5  is  approximately 
bounded  by  the  circular  arc: 

(^2  - 0.1)^  + (7r3  - 0.1)^  = (0.25  - 0.1)^  , 
for  the  case  in  which  n = 30  and  a = 4. 
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TABLE  1.  SUMMARY  STATISTICS  FROM  THE  DISTRIBUTION  OF  THE  DISCRETE  RANGE 
FROM  THREE  SETS  OF  BERNOULLI  TRIALS 
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Figure  3.  Operating  Characteristics  for  Three  Non-Parametric  Tests 
Based  Upon  the  Binomial  Range  of  Three  Sets  of  Bernoulli  Trials 


1 


accept  given  Ti  =0*1, 


Figure  4.  Two-Dimensional  Operating  Characteristic  for  a Non- Parametric  Test  Using  the 
Binomial  Range  from  Three  Sets  of  Bernoulli  Trials  (Sample  Size  30) 


ATTACHMENT  2 


COMPUTER  SOURCE  PROGRAMS  FOR  OBTAINING  THE 

PROBABILITY  DISTRIBUTION  FUNCTION  OF  THE  DISCRETE  BINOMIAL 
RANGE  FOR  THREE  SETS  OF  BERNOULLI  TRIALS 

Three  source  programs  are  given:  an  executive  program  for  I/O 
and  subprogram  calls,  MAIN;  a subroutine  for  computing  the  distribution 
of  binomial  range,  BINRNG;  and  a function  for  calculating  the  binomial 
probability,  PBERN.  All  programs  are  written  in  FORTRAN  4 for  the  IBM 
360  computer. 

Input  requirements  are:  (1)  an  alphameric  title  card  and  (2)  a 
card  specifying  the  sample  size  and  binomial  probability  parameters  for 
each  of  the  three  experiments.  Output  echoes  input  and  lists  the  p.d.f., 
c.d.f.,  and  upper  tail  probabilities  for  the  discrete  range.  An  example 
is  provided. 


f '2- 
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XLCUT  IVE  FOR  PBtRN**<*^*^*******<»<*******»****»*<Kn»»** 
Ma1!4  PROGRAN-  to  develop  a set  of  DliiTRIBUTION  FUNCTIONS 

FOP  THE  Range  of  ouicuhes  erom  three  experiments,  each 

COt;SlqTI.;G  OF  N 'jERMOULLl  iPlALS 
IMPLICIT  REAL^^e  (A-H,0-Z) 

D I PENSION  TirLF(2  0),PS(100)  , CDF  ( 1 ()  0 >.,  PMEAN  ( 3 ) ,N5AMP(3) 

1 COf>,TrJUE 

RFaO  (5»  lOri  ,EmD  = 30)  T ITLE,NSAMP  ( 1 ) ,NSaMP  (2)  ,NSAMP  (3)  , 

1 P^  EaN  ( 1 ) ,Pi'LAN  (2)  ,PEi£AN  (3) 

100  FOPf-^T  (2((A4/3I3.  lX,3F10.n) 

■-•.'EUTt  (6,200)  TirLE,r.SA.'-1P(l)  ,NSAMP(2)  ,NSaMP(3)  , 

1 PMF  AN  ( 1 ) tf’MEAN  (2)  ,Pa>Ean  (3) 

200  FOf'-'AT  ( iHl  ,20AA/1HO,  *SamhLE  SIZES  aRE  ; * , 3 ( 3X  , 1 3 ) , 

1 * M7H  trial  PR0;^S.  ; » , T (2X,F10.E)  ) 

IJ  = :iA/0  (NSAMP  ( 1 ) ,NSAMP  (2)  ,i’J5AMP  (3)  ) 

NP  ] =i\'  4-  1 

WRITE  HcAUlMGS 
i^RITe:  (E),30n) 

300  FOkMaT  ( 1HO,7X,3H:jO.  ,3X,7HL)E  NSITY,  1 IH  CUMUL.  PR.,9H  REM,  PR.) 

Call  p i NK^Jo  ( hs amp  , pme'  an  , ps  , cdf  , i o o ) 

00  7 1=1,NP1 

1 M 1 = i - 1 

RuE  =1  . OuO-CUf  ( I ) 

WPITE  (6,400)  IMl ,PS ( I ) ,CDE  ( I ) ,ROF 
400  format ( /X, 13,3F10,4) 

7 COiMl'lue 
SU:  il=0.  QUO 
SUM2  = 0.()U0 
DO  b I = 1,'1P1 
FlrDFLOA  r ( I-l ) 

SU)‘1=5UM1  ♦Flops  ( I ) 

SUE.2  = SUM'2  + FI<>F  I<>PS  ( I ) 

5 CON  T I /jUE 

VARS  = SUM2-SUlir-><>2 
ST['f.'.V=nbDKT  (VARS) 

COE  V,'.  =bTi)DV/SU'1l 
WRITE  (6,10)  5Ui;l  ,STDDV,COFVN 
10  FOM'IaT  ( IMO  ,9X,6HMEAN.R,  1 ()X,5HSr[)DV  ,t3».  ,7hC0F  VAR/3F15.5) 

GO  TO  1 
30  C'jr;TirJUE 
CALL  EXIT 
STOP 
E.\-0 
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SU'.POJT  lUE  filKRNb{NSAHP,PMEANfPS»CDf"  trJDlM) 

C 

C SUF’ROUT  ir;t  TO  OHTaIU  THE  [)  I S T R I riU  T I Oh..  FUNCTIOIm'  OF  THE  RANGE  OF  A 

0 dir^Of.iAL  vahiamle  From  Three  in[)epfj.n ^ent  experiments*  each 

c OF  WHICH  Co.NSiSTS  OF  N BERNOULLI  TRIALS  WITH  PROB,  PMEaN, 

IMPLICIT  REaL^P  (A“H,0-Z) 

INTEGER  R,5  . ' 

OIN.ENSlOfi  PS(NOIM)  ,CDF{NL»IM)  ,PMEAN(-J)  ,NSAMP(3) 

C 

C NSAMP  SAM^^LE  SIZE  OF  THE  I TH  EXPERIMENTAL  SET 

C N ^vx  sample  size  or  HERUOULLI  EXPE-uhLiNTS 

C PS(f!S)^  PROOAHILirY  DEN.SITY  EONCTlO,,i  oF  THE  RANGE  FROM  THREE  EXP’mTS 

C CDFC'S)  THE  CU'MjLATIVE  l>  I 3TR  I BU  T I Qi'i  FUNCTION  OF  THE  ABOVE  RANGE 

C PoER'' (K  ,:•J,P^■.^A''l)  IS  A FUNCTIOrJ  whICH  CALCULATES  THE  BINOMIAL 

C P.U.K.  H'R  At'blJfiENT  K «’ITH  PARAMcTEN^b  M — THE  SAMPLE  SIZE  — 

C ANO  f'Hf3 Arc-THE  PROBABILITY  OF  THE  E^'ElgT  (SUCCESS)  ON  A SINGLE  TRIAL. 

C K IS  A MEMBLR  OK  Trtt  SET  {0*N). 

C OROE'-  sample  size  FROM  LARIYEST  TO  SMALLEST 

Du  BG  1=1 tP 
IP1=I+1 

DO  3(WI  = 1F1,3 

IF  (.••iSAf.P  (1  ) ...’E.nSa.iP  { I I ) ) GO  TO  31 
NH=NSAnP ( I ) 

MSaMP  ( I ) =NSAriP  (ID 
NSAMf  (ID  =Nii 
HuL('  = RHE  AN  { I ) 

PMEAr:  { I ) =Pt'.E  an  (ID 
Pi-iEAf  (II)=HOLO 
31  COI.rTE'Ut 
30  C(.'OTri;UE 
20  COOT  1 nut 
N = ,nSa''’P(1) 

c 

c INITIALIZE  THE  R ANGE- ARGUME NT  (NS)  LuOP 

C 

TOTAL =0.000 
NPl=f;*  1 

C start  ',S  LOOi^ 

DO  1 iJS=ltNPl 
S=nS-l 
PSUM=0.0D0 
C 

C start  summation  LOOP 

c 

DO  2 NR=1,NP1 
R = Nfv-l 
C 

c CALCULaTL  THE  FIRST  FACTOR 

C 

Fl=Pt-EPI](R,N,PMEAN{l)  ) 


C 
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n n o non 


I 


-c  SLLE.CT  The  Terms  of  the  second  factor 

c 

K1=R-S 
K£'  = R4S 

IF(KI.LT.O)  uO  to  3 
Tl=Pf>e:PN(Kl  *riSAMP(2)  ,PMEAN(2)  ) 

TTl=P(JERri  (K1  ,(;SAMP  (3)  .PMEAN  (3)  ) 

GO  TO  A 
3 Ti=0.U0 
TT 1 =0 . njO 

A IF (K2.GT .N.UR.S.EO. 0)  GO  TO  5 
T2=pf'EPi\  (K2,')SAM1>  (2)  fPf  EAN(2)  ) 

TT2=PE-'EK.!  (K2  ,:1SAMP  <3)  »PMEA;'i(3)  ) 

GO  To  0 
5 T2  = 0.i,'n 
TT2=o,ni'n 
b CONTINUE 
F?  = T] ♦ 12 
Ff 2 = TT1^  rT2 

OEFIUE  ImE  HANGE  limits  OF  THE  THIRD  FACTOR 
LLOv,  = MaXO  ( 0 »K-S) 

LUpK  = '-,IlgO  (^i,R♦S) 

IF  (LL  Oi'..LQ.O.Ar-)L'.LUPR.EU.N)  GO  TO  6 ' 

F3  = 0.L)0 
FF 3=0, 000 
LL0=LL0a4 1 
LUP  = l.UPR*l 
00  7 r=LLO»LUP 
K3=K-1 

F 3 = F3*Pl'EKM  (KG^MSAMp  (3)  » PMEAN  ( 3)  ) 

F F' 3 =F F 3* Pfj LKi  j ( K3  « :J 5 AMP  ( 2 ) » PMEaN  ( 2 ) ) 

7 COriTjrjUE 
GO  To  0 
B F3= 1.000 
F F3=l .000 
9 CH:-T  jtjDL 
F 1=F  1<>1 .020 
F 2=F2oi .020 
F 3=F  3-»l  .020 
F F2=F  F2'M  .o2o 
FF  3=FF3<'1 .1)20 

PSUM  = PSUF  ♦Fl<MF2<»F3tFF2<>FF3-F2i>FT2) 

2 COO'TinuE 

END  OF  SUMN'aIIOn  LOOP;  FILL  THE  PRObABILITY  DENSITY  VECTOR  PS, 

PS (NS) =PSUM«i ,0-60 
total  =TOTaL  + FS (MS) 
iFd.OOO-rOTAL.LE, 1,00-5)  GO  TO  11 
1 CONTINUE 
11  COfiTlNUL 
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rt)  i=\bPi«NPi 
ps(  I ) =n.0f'0 
6 CO.\Tl  U£ 

tlxU  OF  PANGF-ARGUMENT  LOOP;  DEVELOP  ThF;  CUMULATIVE  DISTRIBUTION  FUNCTION 

CCF ( 1 ) =PS  ( 1 ) 

DC  In  i=E«f;pi 

CDF ( 1 ) =CCF ( I-l ) *05 ( I ) 

0 C0\T1.MjL 


o o n n o 


>>>>>>>>  PBERN 


C>>>>>>>>>>>>>> 

— Fu.''CTro';  suukoutkje  to  calculate  the  Value  of  the  binomial  proba- 
pility 

SIlviOMIAL  coefficient  {N»K)  * (P)<>*K  * (1-P)«*(N-K) 

FUNCTION  PULkN  (K»fi»P) 

Implicit  ( A-h»U“Z) 

PtEH-i=n.ouo 

IF(K.GT.N)  f^hTURN 

0 = - P 

PhEP'l  = u<n»r.  , 

IF  (K.KU.O)  RETUFLN 
C.Hj  ir  1 = 1 ^t< 

10  = PoEkU  (UFL0AT{M-I^1)"P)  / (UFLOATd)  « 0) 

RFTUFN 

END 
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OISTKIhUT  lO.'l  oF  The  DISCRETE  FKQM  TmREE  SlTS  OF  rfEWNOULLi  TRIALS 


SIZES 

are: 

AO  30 

30  WI 

NO. 

UEJSITY 

CU-<UL. 

KtM.  PR 

{) 

0.00t>b 

O.OObS 

0.993b 

1 

0.0S62 

0.0629 

0.9374 

2 

0.1233 

0.1960 

0.H140 

3 

0 • 1 76d 

0.3627 

0.6373 

*4 

0. lH6b 

0.4bl)7 

b 

U* Ibdl 

0.7073 

0.2927 

d 

0.1 1d3 

C.' 236 

0.1764 

7 

0.0777 

n . 0 1 3 

0.0967 

3 

0.0477 

o.nsio 

9 

0.020b 

n.'-)756 

0.0242 

U‘ 

0.0137 

0 . 9395 

0.0  lob 

11 

0 . 00b3 

O.'i'JSV 

0.0041 

12 

0.0027 

(1.  19Ab 

0.001b 

13 

0.0010 

0.9995 

0.000b 

U> 

0.0003 

0 .9-v09 

0.0001 

lb 

O.OOul 

1 . nooO 

0 . 0 0 0 0 

16 

0. 0000 

1 . 1'  0 n () 

u . 1)  0 0 0 

1 / 

0.0 

1 . 0 0 n 0 

u . oooo 

IS 

0.0 

1 .onnO 

0,00  no 

19 

0.0 

1 . fi  0 0 0 

C . 0 0 0 0 

20 

0.0 

1 . nniH) 

0 . n 0 n 0 

21 

0.0 

1 . n 0 0 u 

0.  oooo 

22 

0.0 

1 . 1'ono 

0 . 0 0 0 0 

23 

c.o 

1 . n 0 0 0 

0 . oooo 

2^ 

0.0 

1 . n 0 n 0 

0 . 0 (;  0 0 

2b 

0.0 

) . r 0 0 0 

0.0000 

2d 

0.0 

1 . r 0 0 j 

0 . 0 0 0 0 

27 

0.0 

1 .0000 

. 0 0 0 0 

2E 

0 . 0 

1 . rn.oo 

.000  0 

29 

u.o 

1 . D 0 0 0 

. 0 n 0 0 

30 

0.0 

1 . 0 0 n 0 

0 . oooo 

31 

0.0 

1 . oooo 

n.oooo 

32 

0.0 

l . 0000 

0 . 0 0 0 0 

33 

0.0 

1.00  0 0 

0.0000 

3^ 

0.0 

1 . i.OOO 

0.0000 

3b 

0.0 

l . oooo 

f) . 0 0 n 0 

3d 

0.0 

1 . 0 0 0 0 

o.oiuio 

37 

0.0 

1 . 0 0 0 0 

O.OGOO 

3m 

0.0 

1.0000 

0 . i''  0 0 0 

39 

0.0 

1 . 0 0 0 0 

0.0000 

4 0 

0.0 

1 . 0, 0 0 0 

0 . 0(>00 

‘'LAi\,R 

A923H 


STOjV 

2.2129^ 


COF  VAR 
0.^9260 


0.1000  0.1000 


0.2500 
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